12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415 |
- function varargout = spm_sp(varargin)
- % Orthogonal (design) matrix space setting & manipulation
- % FORMAT varargout = spm_spc(action,varargin)
- %
- % This function computes the different projectors related to the row
- % and column spaces X. It should be used to avoid redundant computation
- % of svd on large X matrix. It is divided into actions that set up the
- % space, (Create,Set,...) and actions that compute projections (pinv,
- % pinvXpX, pinvXXp, ...) This is motivated by the problem of rounding
- % errors that can invalidate some computation and is a tool to work
- % with spaces.
- %
- % The only thing that is not easily computed is the null space of
- % the line of X (assuming size(X,1) > size(X,2)).
- % To get this space (a basis of it or a projector on it) use spm_sp on X'.
- %
- % The only restriction on the use of the space structure is when X is
- % so big that you can't fit X and its svd in memory at the same time.
- % Otherwise, the use of spm_sp will generally speed up computations and
- % optimise memory use.
- %
- % Note that since the design matrix is stored in the space structure,
- % there is no need to keep a separate copy of it.
- %
- % ----------------
- %
- % The structure is:
- % x = struct(...
- % 'X', [],... % Mtx
- % 'tol', [],... % tolerance
- % 'ds', [],... % vectors of singular values
- % 'u', [],... % u as in X = u*diag(ds)*v'
- % 'v', [],... % v as in X = u*diag(ds)*v'
- % 'rk', [],... % rank
- % 'oP', [],... % orthogonal projector on X
- % 'oPp', [],... % orthogonal projector on X'
- % 'ups', [],... % space in which this one is embeded
- % 'sus', []); % subspace
- %
- % The basic required fields are X, tol, ds, u, v, rk.
- %
- % ======================================================================
- %
- % FORMAT x = spm_sp('Set',X)
- % Set up space structure, storing matrix, singular values, rank & tolerance
- % X - a (design) matrix (2D)
- % x - the corresponding space structure, with basic fields filled in
- % The SVD is an "economy size" svd, using MatLab's svd(X,0)
- %
- %
- % FORMAT r = spm_sp('oP',x[,Y])
- % FORMAT r = spm_sp('oPp',x[,Y])
- % Return orthogonal projectors, or orthogonal projection of data Y (if passed)
- % x - space structure of matrix X
- % r - ('oP' usage) ortho. projection matrix projecting into column space of x.X
- % - ('oPp' usage) ortho. projection matrix projecting into row space of x.X
- % Y - data (optional)
- % - If data are specified then the corresponding projection of data is
- % returned. This is usually more efficient that computing and applying
- % the projection matrix directly.
- %
- %
- % FORMAT pX = spm_sp('pinv',x)
- % Returns a pseudo-inverse of X - pinv(X) - computed efficiently
- % x - space structure of matrix X
- % pX - pseudo-inverse of X
- % This is the same as MatLab's pinv - the Moore-Penrose pseudoinverse
- % ( Note that because size(pinv(X)) == size(X'), it is not generally )
- % ( useful to compute pinv(X)*Data sequentially (as is the case for )
- % ( 'res' or 'oP') )
- %
- %
- % FORMAT pXpX = spm_sp('pinvxpx',x)
- % Returns a pseudo-inverse of X'X - pinv(X'*X) - computed efficiently
- % x - space structure of matrix X
- % pXpX - pseudo-inverse of (X'X)
- % ( Note that because size(pinv(X'*X)) == [size(X,2) size(X,2)], )
- % ( it is not useful to compute pinv(X'X)*Data sequentially unless )
- % ( size(X,1) < size(X,2) )
- %
- %
- % FORMAT XpX = spm_sp('xpx',x)
- % Returns (X'X) - computed efficiently
- % x - space structure of matrix X
- % XpX - (X'X)
- %
- %
- % FORMAT pXXp = spm_sp('pinvxxp',x)
- % Returns a pseudo-inverse of XX' - pinv(X*X') - computed efficiently
- % x - space structure of matrix X
- % pXXp - pseudo-inverse of (XX')
- %
- %
- % FORMAT XXp = spm_sp('xxp',x)
- % Returns (XX') - computed efficiently
- % x - space structure of matrix X
- % XXp - (XX')
- %
- %
- % FORMAT b = spm_sp('isinsp',x,c[,tol])
- % FORMAT b = spm_sp('isinspp',x,c[,tol])
- % Check whether vectors c are in the column/row space of X
- % x - space structure of matrix X
- % c - vector(s) (Multiple vectors passed as a matrix)
- % tol - (optional) tolerance (for rounding error)
- % [defaults to tolerance specified in space structure: x.tol]
- % b - ('isinsp' usage) true if c is in the column space of X
- % - ('isinspp' usage) true if c is in the column space of X
- %
- % FORMAT b = spm_sp('eachinsp',x,c[,tol])
- % FORMAT b = spm_sp('eachinspp',x,c[,tol])
- % Same as 'isinsp' and 'isinspp' but returns a logical row vector of
- % length size(c,2).
- %
- % FORMAT N = spm_sp('n',x)
- % Simply returns the null space of matrix X (same as matlab NULL)
- % (Null space = vectors associated with zero eigenvalues)
- % x - space structure of matrix X
- % N - null space
- %
- %
- % FORMAT r = spm_sp('nop',x[,Y])
- % Orthogonal projector onto null space of X, or projection of data Y (if passed)
- % x - space structure of matrix X
- % Y - (optional) data
- % r - (if no Y passed) orthogonal projection matrix into the null space of X
- % - (if Y passed ) orthogonal projection of data into the null space of X
- % ( Note that if xp = spm_sp('set',x.X'), we have: )
- % ( spm_sp('nop',x) == spm_sp('res',xp) )
- % ( or, equivalently: )
- % ( spm_sp('nop',x) + spm_sp('oP',xp) == eye(size(xp.X,1)); )
- %
- %
- % FORMAT r = spm_sp('res',x[,Y])
- % Returns residual formaing matrix wirit column space of X, or residuals (if Y)
- % x - space structure of matrix X
- % Y - (optional) data
- % r - (if no Y passed) residual forming matrix for design matrix X
- % - (if Y passed ) residuals, i.e. residual forming matrix times data
- % ( This will be more efficient than
- % ( spm_sp('res',x)*Data, when size(X,1) > size(X,2)
- % Note that this can also be seen as the orthogonal projector onto the
- % null space of x.X' (which is not generally computed in svd, unless
- % size(X,1) < size(X,2)).
- %
- %
- % FORMAT oX = spm_sp('ox', x)
- % FORMAT oXp = spm_sp('oxp',x)
- % Returns an orthonormal basis for X ('ox' usage) or X' ('oxp' usage)
- % x - space structure of matrix X
- % oX - orthonormal basis for X - same as orth(x.X)
- % xOp - *an* orthonormal for X' (but not the same as orth(x.X'))
- %
- %
- % FORMAT b = spm_sp('isspc',x)
- % Check a variable is a structure with the right fields for a space structure
- % x - candidate variable
- % b - true if x is a structure with fieldnames corresponding to spm_sp('create')
- %
- %
- % FORMAT [b,e] = spm_sp('issetspc',x)
- % Test whether a variable is a space structure with the basic fields set
- % x - candidate variable
- % b - true is x is a structure with fieldnames corresponding to
- % spm_sp('Create'), which has it's basic fields filled in.
- % e - string describing why x fails the issetspc test (if it does)
- % This is simply a gateway function combining spm_sp('isspc',x) with
- % the internal subfunction sf_isset, which checks that the basic fields
- % are not empty. See sf_isset (below).
- %
- %-----------------------------------------------------------------------
- % SUBFUNCTIONS:
- %
- % FORMAT b = sf_isset(x)
- % Checks that the basic fields are non-empty (doesn't check they're right!)
- % x - space structure
- % b - true if the basic fields are non-empty
- %_______________________________________________________________________
- % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
- % Jean-Baptiste Poline
- % $Id: spm_sp.m 5219 2013-01-29 17:07:07Z spm $
- if nargin==0
- error('Do what? no arguments given...')
- else
- action = varargin{1};
- end
- %- check the very basics arguments
- switch lower(action),
- case {'create','set','issetspc','isspc'}
- %- do nothing
- otherwise,
- if nargin==1, error('No space : can''t do much!'), end
- [ok,str] = spm_sp('issetspc',varargin{2});
- if ~ok, error(str), else, sX = varargin{2}; end;
- end;
- switch lower(action),
- case 'create' %-Create space structure
- %=======================================================================
- % x = spm_sp('Create')
- varargout = {sf_create};
- case 'set' %-Set singular values, space basis, rank & tolerance
- %=======================================================================
- % x = spm_sp('Set',X)
- if nargin==1, error('No design matrix : can''t do much!'),
- else X = varargin{2}; end
- if isempty(X), varargout = {sf_create}; return, end
- %- only sets plain matrices
- %- check X has 2 dim only
- if max(size(size(X))) > 2, error('Too many dim in the set'), end
- if ~isnumeric(X), error('only sets numeric matrices'), end
- varargout = {sf_set(X)};
- case {'p', 'transp'} %-Transpose space of X
- %=======================================================================
- switch nargin
- case 2
- varargout = {sf_transp(sX)};
- otherwise
- error('too many input argument in spm_sp');
- end % switch nargin
- case {'op', 'op:'} %-Orthogonal projectors on space of X
- %=======================================================================
- % r = spm_sp('oP', sX[,Y])
- % r = spm_sp('oP:', sX[,Y]) %- set to 0 less than tolerence values
- %
- % if isempty(Y) returns as if Y not given
- %-----------------------------------------------------------------------
- switch nargin
- case 2
- switch lower(action),
- case 'op'
- varargout = {sf_op(sX)};
- case 'op:'
- varargout = {sf_tol(sf_op(sX),sX.tol)};
- end %- switch lower(action),
- case 3
- Y = varargin{3};
- if isempty(Y), varargout = {spm_sp(action,sX)}; return, end
- if size(Y,1) ~= sf_s1(sX), error('Dim dont match'); end;
- switch lower(action),
- case 'op'
- varargout = {sf_op(sX)*Y};
- case 'op:'
- varargout = {sf_tol(sf_op(sX)*Y,sX.tol)};
- end % switch lower(action)
- otherwise
- error('too many input argument in spm_sp');
- end % switch nargin
- case {'opp', 'opp:'} %-Orthogonal projectors on space of X'
- %=======================================================================
- % r = spm_sp('oPp',sX[,Y])
- % r = spm_sp('oPp:',sX[,Y]) %- set to 0 less than tolerence values
- %
- % if isempty(Y) returns as if Y not given
- %-----------------------------------------------------------------------
- switch nargin
- case 2
- switch lower(action),
- case 'opp'
- varargout = {sf_opp(sX)};
- case 'opp:'
- varargout = {sf_tol(sf_opp(sX),sX.tol)};
- end %- switch lower(action),
- case 3
- Y = varargin{3};
- if isempty(Y), varargout = {spm_sp(action,sX)}; return, end
- if size(Y,1) ~= sf_s2(sX), error('Dim dont match'); end;
- switch lower(action),
- case 'opp'
- varargout = {sf_opp(sX)*Y};
- case 'opp:'
- varargout = {sf_tol(sf_opp(sX)*Y,sX.tol)};
- end % switch lower(action)
- otherwise
- error('too many input argument in spm_sp');
- end % switch nargin
- case {'x-','x-:'} %-Pseudo-inverse of X - pinv(X)
- %=======================================================================
- % = spm_sp('x-',x)
- switch nargin
- case 2
- switch lower(action),
- case {'x-'}
- varargout = { sf_pinv(sX) };
- case {'x-:'}
- varargout = {sf_tol( sf_pinv(sX), sf_t(sX) )};
- end
- case 3
- %- check dimensions of Y
- Y = varargin{3};
- if isempty(Y), varargout = {spm_sp(action,sX)}; return, end
- if size(Y,1) ~= sf_s1(sX), error(['Dim dont match ' action]); end
- switch lower(action),
- case {'x-'}
- varargout = { sf_pinv(sX)*Y };
- case {'x-:'}
- varargout = {sf_tol( sf_pinv(sX)*Y, sf_t(sX) )};
- end
- otherwise
- error(['too many input argument in spm_sp ' action]);
- end % switch nargin
- case {'xp-','xp-:','x-p','x-p:'} %- Pseudo-inverse of X'
- %=======================================================================
- % pX = spm_sp('xp-',x)
- switch nargin
- case 2
- switch lower(action),
- case {'xp-','x-p'}
- varargout = { sf_pinvxp(sX) };
- case {'xp-:','x-p:'}
- varargout = {sf_tol( sf_pinvxp(sX), sf_t(sX) )};
- end
- case 3
- %- check dimensions of Y
- Y = varargin{3};
- if isempty(Y), varargout = {spm_sp(action,sX)}; return, end
- if size(Y,1) ~= sf_s2(sX), error(['Dim dont match ' action]); end
- switch lower(action),
- case {'xp-','x-p'}
- varargout = { sf_pinvxp(sX)*Y };
- case {'xp-:','x-p:'}
- varargout = {sf_tol( sf_pinvxp(sX)*Y, sf_t(sX) )};
- end
- otherwise
- error(['too many input argument in spm_sp ' action]);
- end % switch nargin
- case {'cukxp-','cukxp-:'} %- Coordinates of pinv(X') in the base of uk
- %=======================================================================
- % pX = spm_sp('cukxp-',x)
- switch nargin
- case 2
- switch lower(action),
- case {'cukxp-'}
- varargout = { sf_cukpinvxp(sX) };
- case {'cukxp-:'}
- varargout = {sf_tol(sf_cukpinvxp(sX),sX.tol)};
- end
- case 3
- %- check dimensions of Y
- Y = varargin{3};
- if isempty(Y), varargout = {spm_sp(action,sX)}; return, end
- if size(Y,1) ~= sf_s2(sX), error(['Dim dont match ' action]); end
- switch lower(action),
- case {'cukxp-'}
- varargout = { sf_cukpinvxp(sX)*Y };
- case {'cukxp-:'}
- varargout = {sf_tol(sf_cukpinvxp(sX)*Y,sX.tol)};
- end
- otherwise
- error(['too many input argument in spm_sp ' action]);
- end % switch nargin
- case {'cukx','cukx:'} %- Coordinates of X in the base of uk
- %=======================================================================
- % pX = spm_sp('cukx',x)
- switch nargin
- case 2
- switch lower(action),
- case {'cukx'}
- varargout = { sf_cukx(sX) };
- case {'cukx:'}
- varargout = {sf_tol(sf_cukx(sX),sX.tol)};
- end
- case 3
- %- check dimensions of Y
- Y = varargin{3};
- if isempty(Y), varargout = {spm_sp(action,sX)}; return, end
- if size(Y,1) ~= sf_s2(sX), error(['Dim dont match ' action]); end
- switch lower(action),
- case {'cukx'}
- varargout = { sf_cukx(sX)*Y };
- case {'cukx:'}
- varargout = {sf_tol(sf_cukx(sX)*Y,sX.tol)};
- end
- otherwise
- error(['too many input argument in spm_sp ' action]);
- end % switch nargin
- case {'rk'} %- Returns rank
- %=======================================================================
- varargout = { sf_rk(sX) };
- case {'ox', 'oxp'} %-Orthonormal basis sets for X / X'
- %=======================================================================
- % oX = spm_sp('ox', x)
- % oXp = spm_sp('oxp',x)
- if sf_rk(sX) > 0
- switch lower(action)
- case 'ox'
- varargout = {sf_uk(sX)};
- case 'oxp'
- varargout = {sf_vk(sX)};
- end
- else
- switch lower(action)
- case 'ox'
- varargout = {zeros(sf_s1(sX),1)};
- case 'oxp'
- varargout = {zeros(sf_s2(sX),1)};
- end
- end
- case {'x', 'xp'} %- X / X' robust to spm_sp changes
- %=======================================================================
- % X = spm_sp('x', x)
- % X' = spm_sp('xp',x)
- switch lower(action)
- case 'x', varargout = {sX.X};
- case 'xp', varargout = {sX.X'};
- end
- case {'xi', 'xpi'} %- X(:,i) / X'(:,i) robust to spm_sp changes
- %=======================================================================
- % X = spm_sp('xi', x)
- % X' = spm_sp('xpi',x)
- i = varargin{3}; % NO CHECKING on i !!! assumes correct
- switch lower(action)
- case 'xi', varargout = {sX.X(:,i)};
- case 'xpi', varargout = {sX.X(i,:)'};
- end
- case {'uk','uk:'} %- Returns u(:,1:r)
- %=======================================================================
- % pX = spm_sp('uk',x)
- % Notice the difference with 'ox' : 'ox' always returns a basis of the
- % proper siwe while this returns empty if rank is null
- warning('can''t you use ox ?');
- switch nargin
- case 2
- switch lower(action),
- case {'uk'}
- varargout = { sf_uk(sX) };
- case {'uk:'}
- varargout = { sf_tol(sf_uk(sX),sX.tol) };
- end
- case 3
- %- check dimensions of Y
- Y = varargin{3};
- if isempty(Y), varargout = {spm_sp(action,sX)}; return, end
- if size(Y,1) ~= sf_rk(sX), error(['Dim dont match ' action]); end
- switch lower(action),
- case {'uk'}
- varargout = { sf_uk(sX)*Y };
- case {'uk:'}
- varargout = {sf_tol(sf_uk(sX)*Y,sX.tol)};
- end
- otherwise
- error(['too many input argument in spm_sp ' action]);
- end % switch nargin
- case {'pinvxpx', 'xpx-', 'pinvxpx:', 'xpx-:',} %- Pseudo-inv of (X'X)
- %=======================================================================
- % pXpX = spm_sp('pinvxpx',x [,Y])
- switch nargin
- case 2
- switch lower(action),
- case {'xpx-','pinvxpx'}
- varargout = {sf_pinvxpx(sX)};
- case {'xpx-:','pinvxpx:'}
- varargout = {sf_tol(sf_pinvxpx(sX),sX.tol)};
- end %-
- case 3
- %- check dimensions of Y
- Y = varargin{3};
- if isempty(Y), varargout = {spm_sp(action,sX)}; return, end
- if size(Y,1) ~= sf_s2(sX), error('Dim dont match'); end;
- switch lower(action),
- case {'xpx-','pinvxpx'}
- varargout = {sf_pinvxpx(sX)*Y};
- case {'xpx-:','pinvxpx:'}
- varargout = {sf_tol(sf_pinvxpx(sX)*Y,sX.tol)};
- end %-
- otherwise
- error('too many input argument in spm_sp');
- end % switch nargin
- case {'xpx','xpx:'} %-Computation of (X'*X)
- %=======================================================================
- % XpX = spm_sp('xpx',x [,Y])
- switch nargin
- case 2
- switch lower(action),
- case {'xpx'}
- varargout = {sf_xpx(sX)};
- case {'xpx:'}
- varargout = {sf_tol(sf_xpx(sX),sX.tol)};
- end %-
- case 3
- %- check dimensions of Y
- Y = varargin{3};
- if isempty(Y), varargout = {spm_sp(action,sX)}; return, end
- if size(Y,1) ~= sf_s2(sX), error('Dim dont match'); end;
- switch lower(action),
- case {'xpx'}
- varargout = {sf_xpx(sX)*Y};
- case {'xpx:'}
- varargout = {sf_tol(sf_xpx(sX)*Y,sX.tol)};
- end %-
- otherwise
- error('too many input argument in spm_sp');
- end % switch nargin
- case {'cx->cu','cx->cu:'} %-coordinates in the basis of X -> basis u
- %=======================================================================
- %
- % returns cu such that sX.X*cx == sX.u*cu
- switch nargin
- case 2
- switch lower(action),
- case {'cx->cu'}
- varargout = {sf_cxtwdcu(sX)};
- case {'cx->cu:'}
- varargout = {sf_tol(sf_cxtwdcu(sX),sX.tol)};
- end %-
- case 3
- %- check dimensions of Y
- Y = varargin{3};
- if isempty(Y), varargout = {spm_sp(action,sX)}; return, end
- if size(Y,1) ~= sf_s2(sX), error('Dim dont match'); end;
- switch lower(action),
- case {'cx->cu'}
- varargout = {sf_cxtwdcu(sX)*Y};
- case {'cx->cu:'}
- varargout = {sf_tol(sf_cxtwdcu(sX)*Y,sX.tol)};
- end %-
- otherwise
- error('too many input argument in spm_sp');
- end % switch nargin
- case {'xxp-','xxp-:','pinvxxp','pinvxxp:'} %-Pseudo-inverse of (XX')
- %=======================================================================
- % pXXp = spm_sp('pinvxxp',x [,Y])
- switch nargin
- case 2
- switch lower(action),
- case {'xxp-','pinvxxp'}
- varargout = {sf_pinvxxp(sX)};
- case {'xxp-:','pinvxxp:'}
- varargout = {sf_tol(sf_pinvxxp(sX),sX.tol)};
- end %-
- case 3
- %- check dimensions of Y
- Y = varargin{3};
- if isempty(Y), varargout = {spm_sp(action,sX)}; return, end
- if size(Y,1) ~= sf_s1(sX), error('Dim dont match'); end;
- switch lower(action),
- case {'xxp-','pinvxxp'}
- varargout = {sf_pinvxxp(sX)*Y};
- case {'xxp-:','pinvxxp:'}
- varargout = {sf_tol(sf_pinvxxp(sX)*Y,sX.tol)};
- end %-
- otherwise
- error('too many input argument in spm_sp');
- end % switch nargin
-
- case {'xxp','xxp:'} %-Computation of (X*X')
- %=======================================================================
- % XXp = spm_sp('xxp',x)
- switch nargin
- case 2
- switch lower(action),
- case {'xxp'}
- varargout = {sf_xxp(sX)};
- case {'xxp:'}
- varargout = {sf_tol(sf_xxp(sX),sX.tol)};
- end %-
- case 3
- %- check dimensions of Y
- Y = varargin{3};
- if isempty(Y), varargout = {spm_sp(action,sX)}; return, end
- if size(Y,1) ~= sf_s1(sX), error('Dim dont match'); end;
- switch lower(action),
- case {'xxp'}
- varargout = {sf_xxpY(sX,Y)};
- case {'xxp:'}
- varargout = {sf_tol(sf_xxpY(sX,Y),sX.tol)};
- end %-
- otherwise
- error('too many input argument in spm_sp');
- end % switch nargin
- case {'^p','^p:'} %-Computation of v*(diag(s.^n))*v'
- %=======================================================================
- switch nargin
- case {2,3}
- if nargin==2, n = 1; else n = varargin{3}; end;
- if ~isnumeric(n), error('~isnumeric(n)'), end;
- switch lower(action),
- case {'^p'}
- varargout = {sf_jbp(sX,n)};
- case {'^p:'}
- varargout = {sf_tol(sf_jbp(sX,n),sX.tol)};
- end %-
- case 4
- n = varargin{3};
- if ~isnumeric(n), error('~isnumeric(n)'), end;
- Y = varargin{4};
- if isempty(Y), varargout = {spm_sp(action,sX,n)}; return, end
- if size(Y,1) ~= sf_s2(sX), error('Dim dont match'); end;
- switch lower(action),
- case {'^p'}
- varargout = {sf_jbp(sX,n)*Y};
- case {'^p:'}
- varargout = {sf_tol(sf_jbp(sX,n)*Y,sX.tol)};
- end %-
- otherwise
- error('too many input argument in spm_sp');
- end % switch nargin
- case {'^','^:'} %-Computation of v*(diag(s.^n))*v'
- %=======================================================================
- switch nargin
- case {2,3}
- if nargin==2, n = 1; else n = varargin{3}; end;
- if ~isnumeric(n), error('~isnumeric(n)'), end;
- switch lower(action),
- case {'^'}
- varargout = {sf_jb(sX,n)};
- case {'^:'}
- varargout = {sf_tol(sf_jb(sX,n),sX.tol)};
- end %-
- case 4
- n = varargin{3};
- if ~isnumeric(n), error('~isnumeric(n)'), end;
- Y = varargin{4};
- if isempty(Y), varargout = {spm_sp(action,sX,n)}; return, end
- if size(Y,1) ~= sf_s1(sX), error('Dim dont match'); end;
- switch lower(action),
- case {'^'}
- varargout = {sf_jbY(sX,n,Y)};
- case {'^:'}
- varargout = {sf_tol(sf_jbY(sX,n,Y),sX.tol)};
- end %-
- otherwise
- error('too many input argument in spm_sp');
- end % switch nargin
- case {'n'} %-Null space of sX
- %=======================================================================
- switch nargin
- case 2
- varargout = {sf_n(sX)};
- otherwise
- error('too many input argument in spm_sp');
- end % switch nargin
- case {'np'} %-Null space of sX'
- %=======================================================================
- switch nargin
- case 2
- varargout = {sf_n(sf_transp(sX))};
- otherwise
- error('too many input argument in spm_sp');
- end % switch nargin
- case {'nop', 'nop:'} %- project(or)(ion) into null space
- %=======================================================================
- %
- %
- %
- switch nargin
- case 2
- switch lower(action),
- case {'nop'}
- n = sf_n(sX);
- varargout = {n*n'};
- case {'nop:'}
- n = sf_n(sX);
- varargout = {sf_tol(n*n',sX.tol)};
- end %-
- case 3
- %- check dimensions of Y
- Y = varargin{3};
- if isempty(Y), varargout = {spm_sp(action,sX)}; return, end
- if size(Y,1) ~= sf_s2(sX), error('Dim dont match'); end;
- switch lower(action),
- case {'nop'}
- n = sf_n(sX);
- varargout = {n*(n'*Y)};
- case {'nop:'}
- n = sf_n(sX);
- varargout = {sf_tol(n*(n'*Y),sX.tol)};
- end %-
- otherwise
- error('too many input argument in spm_sp');
- end % switch nargin
- case {'nopp', 'nopp:'} %- projector(ion) into null space of X'
- %=======================================================================
- %
- %
- switch nargin
- case 2
- switch lower(action),
- case {'nopp'}
- varargout = {spm_sp('nop',sf_transp(sX))};
- case {'nopp:'}
- varargout = {spm_sp('nop:',sf_transp(sX))};
- end %-
- case 3
- switch lower(action),
- case {'nopp'}
- varargout = {spm_sp('nop',sf_transp(sX),varargin{3})};
- case {'nopp:'}
- varargout = {spm_sp('nop:',sf_transp(sX),varargin{3})};
- end %-
- otherwise
- error('too many input argument in spm_sp');
- end % switch nargin
- case {'res', 'r','r:'} %-Residual formaing matrix / residuals
- %=======================================================================
- % r = spm_sp('res',sX[,Y])
- %
- %- 'res' will become obsolete : use 'r' or 'r:' instead
- %- At some stage, should be merged with 'nop'
- switch nargin
- case 2
- switch lower(action)
- case {'r','res'}
- varargout = {sf_r(sX)};
- case {'r:','res:'}
- varargout = {sf_tol(sf_r(sX),sX.tol)};
- end %-
- case 3
- %- check dimensions of Y
- Y = varargin{3};
- if isempty(Y), varargout = {spm_sp(action,sX)}; return, end
- if size(Y,1) ~= sf_s1(sX), error('Dim dont match'); end;
-
- switch lower(action)
- case {'r','res'}
- varargout = {sf_rY(sX,Y)};
- case {'r:','res:'}
- varargout = {sf_tol(sf_rY(sX,Y),sX.tol)};
- end %-
- otherwise
- error('too many input argument in spm_sp');
- end % switch nargin
- case {':'}
- %=======================================================================
- % spm_sp(':',sX [,Y [,tol]])
- %- Sets Y and tol according to arguments
- if nargin > 4
- error('too many input argument in spm_sp');
- else
- if nargin > 3
- if isnumeric(varargin{4}), tol = varargin{4};
- else error('tol must be numeric');
- end
- else
- tol = sX.tol;
- end
- if nargin > 2
- Y = varargin{3}; %- if isempty, returns empty, end
- else
- Y = sX.X;
- end
- end
- varargout = {sf_tol(Y,tol)};
- case {'isinsp', 'isinspp'} %- is in space or is in dual space
- %=======================================================================
- % b = spm_sp('isinsp',x,c[,tol])
- % b = spm_sp('isinspp',x,c[,tol])
- %-Check whether vectors are in row/column space of X
- %-Check arguments
- %-----------------------------------------------------------------------
- if nargin<3, error('insufficient arguments - action,x,c required'), end
- c = varargin{3}; %- if isempty(c), dim wont match exept for empty sp.
- if nargin<4, tol=sX.tol; else, tol = varargin{4}; end
- %-Compute according to case
- %-----------------------------------------------------------------------
- switch lower(action)
- case 'isinsp'
- %-Check dimensions
- if size(sX.X,1) ~= size(c,1)
- warning('Vector dim don''t match col. dim : not in space !');
- varargout = { 0 }; return;
- end
- varargout = {all(all( abs(sf_op(sX)*c - c) <= tol ))};
- case 'isinspp'
- %- check dimensions
- if size(sX.X,2) ~= size(c,1)
- warning('Vector dim don''t match row dim : not in space !');
- varargout = { 0 }; return;
- end
- varargout = {all(all( abs(sf_opp(sX)*c - c) <= tol ))};
- end
- case {'eachinsp', 'eachinspp'} %- each column of c in space or in dual space
- %=======================================================================
- % b = spm_sp('eachinsp',x,c[,tol])
- % b = spm_sp('eachinspp',x,c[,tol])
- %-Check whether vectors are in row/column space of X
- %-Check arguments
- %-----------------------------------------------------------------------
- if nargin<3, error('insufficient arguments - action,x,c required'), end
- c = varargin{3}; %- if isempty(c), dim wont match exept for empty sp.
- if nargin<4, tol=sX.tol; else, tol = varargin{4}; end
- %-Compute according to case
- %-----------------------------------------------------------------------
- switch lower(action)
- case 'eachinsp'
- %-Check dimensions
- if size(sX.X,1) ~= size(c,1)
- warning('Vector dim don''t match col. dim : not in space !');
- varargout = { 0 }; return;
- end
- varargout = {all( abs(sf_op(sX)*c - c) <= tol )};
- case 'eachinspp'
- %- check dimensions
- if size(sX.X,2) ~= size(c,1)
- warning('Vector dim don''t match row dim : not in space !');
- varargout = { 0 }; return;
- end
- varargout = {all( abs(sf_opp(sX)*c - c) <= tol )};
- end
- case '==' % test wether two spaces are the same
- %=======================================================================
- % b = spm_sp('==',x1,X2)
- if nargin~=3, error('too few/many input arguments - need 3');
- else X2 = varargin{3}; end;
- if isempty(sX.X)
- if isempty(X2),
- warning('Both spaces empty');
- varargout = { 1 };
- else
- warning('one space empty');
- varargout = { 0 };
- end;
- else
- x2 = spm_sp('Set',X2);
- maxtol = max(sX.tol,x2.tol);
- varargout = { all( spm_sp('isinsp',sX,X2,maxtol)) & ...
- all( spm_sp('isinsp',x2,sX.X,maxtol) ) };
- %- I have encountered one case where the max of tol was needed.
- end;
- case 'isspc' %-Space structure check
- %=======================================================================
- % [b,str] = spm_sp('isspc',x)
- if nargin~=2, error('too few/many input arguments - need 2'), end
- %-Check we've been passed a structure
- if ~isstruct(varargin{2}), varargout={0}; return, end
- %-Go through required field names checking their existance
- % (Get fieldnames once and compare: isfield doesn't work for multiple )
- % (fields, and repeated calls to isfield would result in unnecessary )
- % (replicate calls to fieldnames(varargin{2}). )
- b = 1;
- fnames = fieldnames(varargin{2});
- for str = fieldnames(sf_create)'
- b = b & any(strcmp(str,fnames));
- if ~b, break, end
- end
- if nargout > 1,
- if b, str = 'ok'; else, str = 'not a space'; end;
- varargout = {b,str};
- else, varargout = {b}; end;
- case 'issetspc' %-Is this a completed space structure?
- %=======================================================================
- % [b,e] = spm_sp('issetspc',x)
- if nargin~=2, error('too few/many input arguments - need 2'), end
- if ~spm_sp('isspc',varargin{2})
- varargout = {0,'not a space structure (wrong fieldnames)'};
- elseif ~sf_isset(varargin{2})
- %-Basic fields aren't filled in
- varargout = {0,'space not defined (use ''set'')'};
- else
- varargout = {1,'OK!'};
- end
- case 'size' %- gives the size of sX
- %=======================================================================
- % size = spm_sp('size',x,dim)
- %
- if nargin > 3, error('too many input arguments'), end
- if nargin == 2, dim = []; else dim = varargin{3}; end
- if ~isempty(dim)
- switch dim
- case 1, varargout = { sf_s1(sX) };
- case 2, varargout = { sf_s2(sX) };
- otherwise, error(['unknown dimension in ' action]);
- end
- else %- assumes want both dim
- switch nargout
- case {0,1}
- varargout = { sf_si(sX) };
- case 2
- varargout = { sf_s1(sX), sf_s2(sX) };
- otherwise
- error(['too many output arg in ' mfilename ' ' action]);
- end
- end
- otherwise
- %=======================================================================
- error(['Invalid action (',action,')'])
- %=======================================================================
- end % (case lower(action))
- %=======================================================================
- %- S U B - F U N C T I O N S
- %=======================================================================
- %
- % The rule I tried to follow is that the space structure is accessed
- % only in this sub function part : any sX.whatever should be
- % prohibited elsewhere .... still a lot to clean !!!
- function x = sf_create
- %=======================================================================
- x = struct(...
- 'X', [],... % Matrix
- 'tol', [],... % tolerance
- 'ds', [],... % vectors of singular values
- 'u', [],... % u as in X = u*diag(ds)*v'
- 'v', [], ... % v as in X = u*diag(ds)*v'
- 'rk', [],... % rank
- 'oP', [],... % orthogonal projector on X
- 'oPp', [],... % orthogonal projector on X'
- 'ups', [],... % space in which this one is embeded
- 'sus', []); % subspace
- function x = sf_set(X)
- %=======================================================================
- x = sf_create;
- x.X = X;
- %-- Compute the svd with svd(X,0) : find all the singular values of x.X
- %-- SVD(FULL(A)) will usually perform better than SVDS(A,MIN(SIZE(A)))
- %- if more column that lines, performs on X'
- if size(X,1) < size(X,2)
- [x.v, s, x.u] = svd(full(X'),0);
- else
- [x.u, s, x.v] = svd(full(X),0);
- end
- x.ds = diag(s); clear s;
- %-- compute the tolerance
- x.tol = max(size(x.X))*max(abs(x.ds))*eps;
- %-- compute the rank
- x.rk = sum(x.ds > x.tol);
- function x = sf_transp(x)
- %=======================================================================
- %
- %- Tranpspose the space : note that tmp is not touched, therefore
- %- only contains the address, no duplication of data is performed.
- x.X = x.X';
- tmp = x.v;
- x.v = x.u;
- x.u = tmp;
- tmp = x.oP;
- x.oP = x.oPp;
- x.oPp = tmp;
- clear tmp;
- function b = sf_isset(x)
- %=======================================================================
- b = ~( isempty(x.X) |...
- isempty(x.u) |...
- isempty(x.v) |...
- isempty(x.ds) |...
- isempty(x.tol) |...
- isempty(x.rk) );
- function s1 = sf_s1(x)
- %=======================================================================
- s1 = size(x.X,1);
- function s2 = sf_s2(x)
- %=======================================================================
- s2 = size(x.X,2);
- function si = sf_si(x)
- %=======================================================================
- si = size(x.X);
- function r = sf_rk(x)
- %=======================================================================
- r = x.rk;
- function uk = sf_uk(x)
- %=======================================================================
- uk = x.u(:,1:sf_rk(x));
- function vk = sf_vk(x)
- %=======================================================================
- vk = x.v(:,1:sf_rk(x));
- function sk = sf_sk(x)
- %=======================================================================
- sk = x.ds(1:sf_rk(x));
- function t = sf_t(x)
- %=======================================================================
- t = x.tol;
- function x = sf_tol(x,t)
- %=======================================================================
- x(abs(x) < t) = 0;
- function op = sf_op(sX)
- %=======================================================================
- if sX.rk > 0
- op = sX.u(:,[1:sX.rk])*sX.u(:,[1:sX.rk])';
- else
- op = zeros( size(sX.X,1) );
- end;
- %!!!! to implement : a clever version of sf_opY (see sf_rY)
- function opp = sf_opp(sX)
- %=======================================================================
- if sX.rk > 0
- opp = sX.v(:,[1:sX.rk])*sX.v(:,[1:sX.rk])';
- else
- opp = zeros( size(sX.X,2) );
- end;
- %!!!! to implement : a clever version of sf_oppY (see sf_rY)
- function px = sf_pinv(sX)
- %=======================================================================
- r = sX.rk;
- if r > 0
- px = sX.v(:,1:r)*diag( ones(r,1)./sX.ds(1:r) )*sX.u(:,1:r)';
- else
- px = zeros(size(sX.X,2),size(sX.X,1));
- end
- function px = sf_pinvxp(sX)
- %=======================================================================
- r = sX.rk;
- if r > 0
- px = sX.u(:,1:r)*diag( ones(r,1)./sX.ds(1:r) )*sX.v(:,1:r)';
- else
- px = zeros(size(sX.X));
- end
- function px = sf_pinvxpx(sX)
- %=======================================================================
- r = sX.rk;
- if r > 0
- px = sX.v(:,1:r)*diag( sX.ds(1:r).^(-2) )*sX.v(:,1:r)';
- else
- px = zeros(size(sX.X,2));
- end
- function px = sf_jbp(sX,n)
- %=======================================================================
- r = sX.rk;
- if r > 0
- px = sX.v(:,1:r)*diag( sX.ds(1:r).^(n) )*sX.v(:,1:r)';
- else
- px = zeros(size(sX.X,2));
- end
- function x = sf_jb(sX,n)
- %=======================================================================
- r = sX.rk;
- if r > 0
- x = sX.u(:,1:r)*diag( sX.ds(1:r).^(n) )*sX.u(:,1:r)';
- else
- x = zeros(size(sX.X,1));
- end
- function y = sf_jbY(sX,n,Y)
- %=======================================================================
- r = sX.rk;
- if r > 0
- y = ( sX.u(:,1:r)*diag(sX.ds(1:r).^n) )*(sX.u(:,1:r)'*Y);
- else
- y = zeros(size(sX.X,1),size(Y,2));
- end
- %!!!! to implement : a clever version of sf_jbY (see sf_rY)
- function x = sf_cxtwdcu(sX)
- %=======================================================================
- %- coordinates in sX.X -> coordinates in sX.u(:,1:rk)
- x = diag(sX.ds)*sX.v';
- function x = sf_cukpinvxp(sX)
- %=======================================================================
- %- coordinates of pinv(sX.X') in the basis sX.u(:,1:rk)
- r = sX.rk;
- if r > 0
- x = diag( ones(r,1)./sX.ds(1:r) )*sX.v(:,1:r)';
- else
- x = zeros( size(sX.X,2) );
- end
- function x = sf_cukx(sX)
- %=======================================================================
- %- coordinates of sX.X in the basis sX.u(:,1:rk)
- r = sX.rk;
- if r > 0
- x = diag( sX.ds(1:r) )*sX.v(:,1:r)';
- else
- x = zeros( size(sX.X,2) );
- end
- function x = sf_xpx(sX)
- %=======================================================================
- r = sX.rk;
- if r > 0
- x = sX.v(:,1:r)*diag( sX.ds(1:r).^2 )*sX.v(:,1:r)';
- else
- x = zeros(size(sX.X,2));
- end
- function x = sf_xxp(sX)
- %=======================================================================
- r = sX.rk;
- if r > 0
- x = sX.u(:,1:r)*diag( sX.ds(1:r).^2 )*sX.u(:,1:r)';
- else
- x = zeros(size(sX.X,1));
- end
- function x = sf_xxpY(sX,Y)
- %=======================================================================
- r = sX.rk;
- if r > 0
- x = sX.u(:,1:r)*diag( sX.ds(1:r).^2 )*(sX.u(:,1:r)'*Y);
- else
- x = zeros(size(sX.X,1),size(Y,2));
- end
- function x = sf_pinvxxp(sX)
- %=======================================================================
- r = sX.rk;
- if r > 0
- x = sX.u(:,1:r)*diag( sX.ds(1:r).^(-2) )*sX.u(:,1:r)';
- else
- x = zeros(size(sX.X,1));
- end
- function n = sf_n(sX)
- %=======================================================================
- % if the null space is in sX.v, returns it
- % otherwise, performs Gramm Schmidt orthogonalisation.
- %
- %
- r = sX.rk;
- [q,p]= size(sX.X);
- if r > 0
- if q >= p %- the null space is entirely in v
- if r == p, n = zeros(p,1); else n = sX.v(:,r+1:p); end
- else %- only part of n is in v: same as computing the null sp of sX'
- n = null(sX.X);
- %----- BUG !!!! in that part ----------------------------------------
- %- v = zeros(p,p-q); j = 1; i = 1; z = zeros(p,1);
- %- while i <= p
- %- o = z; o(i) = 1; vpoi = [sX.v(i,:) v(i,1:j-1)]';
- %- o = sf_tol(o - [sX.v v(:,1:j-1)]*vpoi,sX.tol);
- %- if any(o), v(:,j) = o/((o'*o)^(1/2)); j = j + 1; end;
- %- i = i + 1; %- if i>p, error('gramm shmidt error'); end;
- %- end
- %- n = [sX.v(:,r+1:q) v];
- %--------------------------------------------------------------------
- end
- else
- n = eye(p);
- end
- function r = sf_r(sX)
- %=======================================================================
- %-
- %- returns the residual forming matrix for the space sX
- %- for internal use. doesn't Check whether oP exist.
- r = eye(size(sX.X,1)) - sf_op(sX) ;
- function Y = sf_rY(sX,Y)
- %=======================================================================
- % r = spm_sp('r',sX[,Y])
- %
- %- tries to minimise the computation by looking whether we should do
- %- I - u*(u'*Y) or n*(n'*Y) as in 'nop'
- r = sX.rk;
- [q,p]= size(sX.X);
- if r > 0 %- else returns the input;
-
- if r < q-r %- we better do I - u*u'
- Y = Y - sX.u(:,[1:r])*(sX.u(:,[1:r])'*Y); % warning('route1');
- else
- %- is it worth computing the n ortho basis ?
- if size(Y,2) < 5*q
- Y = sf_r(sX)*Y; % warning('route2');
- else
- n = sf_n(sf_transp(sX)); % warning('route3');
- Y = n*(n'*Y);
- end
- end
- end
|