Parafac and Als: Difference between pages

From Eigenvector Research Documentation Wiki
(Difference between pages)
Jump to navigation Jump to search
imported>Jeremy
(Importing text file)
 
imported>Neal
 
Line 1: Line 1:
===Purpose===
===Purpose===


PARAFAC (PARAllel FACtor analysis) for multi-way arrays
Alternating Least Squares computational engine for multivariate curve resolution (MCR).


===Synopsis===
===Synopsis===


:model  = parafac(X,''initval,options'')
:[c,s] = als(x,c0,options);
 
:pred    = parafac(Xnew,model)
 
:options = parafac('options')


===Description===
===Description===


PARAFAC will decompose an array of order ''N'' (where ''N'' ? 3) into the summation over the outer product of ''N'' vectors (a low-rank model). E.g. if ''N''=3 then the array is size ''I'' by ''J'' by ''K''. An example of three-way fluorescence data is shown below..
ALS decomposes a matrix '''X''' as '''CS''' such that '''X''' = '''CS''' + '''E''' where '''E''' is minimized in a least squares sense.  
 
For example, twenty-seven samples containing different amounts of dissolved hydroquinone, tryptophan, phenylalanine, and dopa are measured spectrofluoremetrically using 233 emission wavelengths (250-482 nm) and 24 excitation wavelengths (200-315 nm each 5 nm). A typical sample is also shown.
 
 
A four-component PARAFAC model of these data will give four factors, each corresponding to one of the chemical analytes. This is illustrated graphically below. The first mode scores (loadings in mode 1) in the matrix '''A''' (27�4) contain estimated relative concentrations of the four analytes in the 27 samples. The second mode loadings '''B''' (233�4) are estimated emission loadings and the third mode loadings '''C''' (24�4) are estimated excitation loadings.
 
 
In the PARAFAC algorithm, any missing values must be set to NaN or Inf and are then automatically handled by expectation maximization. This routine employs an alternating least squares (ALS) algorithm in combination with a line search. For 3-way data, the initial estimate of the loadings is usually obtained from the tri-linear decomposition (TLD).
 
====INPUTS====
 
* '''x''' = the multiway array to be decomposed, and
 
* '''ncomp''' =  the number of factors (components) to use, or
 
* '''model''' a PARAFAC model structure (new data are fit to the model i.e. sample mode scores are calculated).
 
====OPTIONAL INPUTS====
 
* '''''initval''''' =  cell array of initial values (initial guess) for the loadings (e.g. model.loads from a previous fit). If not used it can be 0 or [], and
 
*'''''''' options'' =  discussed below.
 
====OUTPUTS====
 
The output model is a structure array with the following fields:
 
* '''modeltype''': 'PARAFAC',
 
* '''datasource''': structure array with information about input data,
 
* '''date''': date of creation,
 
* '''time''': time of creation,
 
* '''info''': additional model information,


* '''loads''': 1 by ''K'' cell array with model loadings for each mode/dimension,
Inputs are the matrix to be decomposed (x) (''M'' by ''N''), and the initial guess (c0). If (c0) is size ''M'' by ''K'', where ''K'' is the number of factors, then it is assumed to be the initial guess for '''C'''. If (c0) is size ''K'' by ''N'' then it is assumed to be the initial guess for '''S'''. If ''M'' = ''K'' then, (c0) is assumed to be the initial guess for '''C'''.


* '''pred''': cell array with model predictions for each input data block,
An optional input (options) is described below.


* '''tsqs''': cell array with T<sup>2</sup> values for each mode,
The outputs are the estimated matrices '''C''', (c), (''M'' by ''K'') and '''S''', (s), (''K'' by ''N''). For a typical decomposition, (c) is a matrix of contributions and (s) is a matrix of spectra. The function


* '''ssqresiduals''': cell array with sum of squares residuals for each mode,
:[c,s] = als(x,c0)


* '''description''': cell array with text description of model, and
will decompose (x) using an non-negatively constrained alternating least squares. To include other constraints, use the (options) input described below.


* '''detail''': sub-structure with additional model details and results.
Note that if no non-zero equality constraints are imposed on a factor the spectra are normalized to unit length. This can lead to significant scaling differences between factors that have non-zero equality constraints and those that do not.
 
The output pred is a structure array that contains the approximation of the data if the options field blockdetails is set to 'all' (see next).


===Options===
===Options===


* '''''options''''' =  a structure array with the following fields:
*     '''display''':  [ 'off' | {'on'} ]  governs level of display to command window.


* '''display''': [ {'on'} | 'off' ], governs level of display,
* '''plots''': [ 'none' | {'final'} ] governs level of plotting.


* '''plots''': [ {'final'} | 'all' | 'none' ], governs level of plotting,
*         '''ccon''': [ 'none' | 'reset' | {'fastnnls'} ] non-negativity on contributions.  (fastnnls = true least-squares solution)


* '''weights''': [], used for fitting a weighted loss function (discussed below),
*         '''scon''': [ 'none' | 'reset' | {'fastnnls'} ] non-negativity on spectra.    (fastnnls = true least-squares solution)


* '''stopcrit''': [1e-6 1e-6 10000 3600] defines the stopping criteria as [(relative tolerance) (absolute tolerance) (maximum number of iterations) (maximum time in seconds)],
*           '''cc''': [ ] contributions equality constraints, must be a matrix with ''M'' rows and up to ''K'' columns with NaN where equality constraints are not applied and real value of the constraint where they are applied. If fewer than ''K'' columns are supplied, the missing columns will be filled in as unconstrained,


* '''init''': [ 0 ], defines how parameters are initialized (discussed below),
*     '''ccwts''': [inf] a scalar value or a 1x''K'' vector with elements corresponding to weightings on constraints (0, no constraint, 0<wt<inf imposes constraint "softly", and inf is hard constrained). If a scalar value is passed for (ccwts), that value is applied for all ''K'' factors.
:Soft constraints are imposed by augmenting an additional column onto '''X''' and a weight of 1 weights a given equality constraint so that it has an equal influence as the average single variable in '''X'''. For additional information see P.J. Gemperline and E. Cash, “Advantages of Soft versus Hard Constraints in Self–Modeling Curve Resolution Problems. Alternating Least Squares with Penalty Functions”, ''Anal. Chem.,'' '''75'''(16), 4236–4243 (2003).


* '''line''': [ 0 | {1}] defines whether to use the line search {default uses it},
*           '''sc''': [ ] spectra equality constraints, must be a matrix with ''N'' columns and up to ''K'' rows with NaN where equality contraints are not applied and real value of the constraint where they are applied.  If fewer than ''K'' rows are supplied, the missing rows will be filled in as unconstrained.


* '''algo''': [ {'ALS'} | 'tld' | 'swatld' ] governs algorithm used,
*       '''scwts''': [inf] weighting for spectral equality constraints (see ccwts) above.


* '''iterative''': settings for iterative reweighted least squares fitting (see help on weights below),
*         '''sclc''': [ ] contributions scale axis, vector with ''M'' elements otherwise 1:M is used.


* '''blockdetails''': 'standard'
*         '''scls''': [ ]  spectra scale axis, vector with ''N'' elements otherwise 1:N is used.


* '''missdat''': this option is not yet active,
*         '''condition''': [{'none'}| 'norm' ] type of conditioning to perform on '''S''' and '''C''' before each regression step. 'norm' uses conditioning which can help stabilize the regression when factors are significantly different in magnitude.


* '''samplemode''': [1], defines which mode should be considered the sample or object mode,
*         '''tolc''': [ {1e-5} ] tolerance on non-negativity for contributions.


* '''constraints''': {3x1 cell}, defines constraints on parameters (discussed below), and
*         '''tols''': [ {1e-5} ]  tolerance on non-negativity for spectra.


* '''coreconsist''': [ {'on'} | 'off' ], governs calculation of core consistency (turning off may save time with large data sets and many components).
*       '''ittol''': [ {1e-8} ] convergence tolerance.


The default options can be retrieved using: options = parafac('options');.
*        '''itmax''':  [ {100} ]  maximum number of iterations.


WEIGHTS
*      '''timemax''':  [ {3600} ]  maximum time for iterations.


Through the use of the ''options'' field weights it is possible to fit a PARAFAC model in a weighted least squares sense The input is an array of the same size as the input data X holding individual weights for each element. The PARAFAC model is then fit in a weighted least squares sense. Instead of minimizing the frobenius norm ||x-M||<sup>2</sup> where M is the PARAFAC model, the norm ||(x-M).\*weights||<sup>2</sup> is minimized. The algorithm used for weighted regression is based on a majorization step according to Kiers, ''Psychometrika'', '''62''', 251-266, 1997 which has the advantage of being computationally inexpensive. If alternatively, the field weights is set to 'iterative' then iteratively reweighted least squares fitting is used. The settings of this can be modified in the field iterative.cutoff_residuals which defines the cutoff for large residuals in terms of the number of robust standard deviations. The lower the number, the more subtle outliers will be ignored.
*    '''rankfail''':  [ 'drop' |{'reset'}| 'random' | 'fail' ]  how are rank deficiencies handled:


INIT
*                  '''drop'''  - drop deficient components from model,


The ''options'' field init is used to govern how the initial guess for the loadings is obtained. If optional input ''initval'' is input then options.init is not used. The following choices for init are available.
*                  '''reset'''  - reset deficient components to initial guess,


Generally, options.init = 0, will do for well-behaved data whereas options.init = 10, will be suitable for difficult models. Difficult models are typically those with many components, with very correlated loadings, or models where there are indications that local minima are present.
*                '''random''' - replace deficient components with random vector,


* '''init''' =  0, PARAFAC chooses initialization {default},
*                 '''fail'''   - stop analysis, give error.


* '''init''' = 1, uses TLD (unless there are missing values then random is used),
===Examples===
 
* '''init''' = 2, initializes loadings with random values,
 
* '''init''' = 3, based on orthogonalization of random values (preferred over 2),
 
* '''init''' = 4, based on singular value decomposition,
 
* '''init''' = 5, based on compression which may be useful for large data, and
 
* '''init''' >  5, based on best fit of many (the value options.init) small runs.
 
CONSTRAINTS
 
The ''options'' field constraints is used to employ constraints on the parameters. It is a cell array with number of elements equal to the number of modes of the input data X. Each cell contains a structure array with the following fields:
 
* '''nonnegativity''': [ {0} | 1 ], a 1 imposes non-negativity.
 
* '''unimodality''': [ {0} | 1 ], a 1 imposes unimodality (1 local maxima).
 
* '''orthogonal''': [ {0} | 1 ], constrain factors in this mode to be orthogonal.
 
* '''orthonormal''': [ {0} | 1 ], constrain factors in this mode to be orthonormal.
 
* '''exponential''': [ {0} | 1 ], a 1 fits an exponential function to the factors in this mode.
 
:smoothness.weight: [0 to 1], imposes smoothness using B-splines, values near 1 impose high smoothness and values close to 0, impose less smoothness.
 
* '''fixed.position''': [ ], a matrix containing 1's and 0's of the same size as the corresponding loading matrix, with a 1 indicating where parameters are fixed.
 
* '''fixed.value''': [ ], a vector containing the fixed values. Thus, if B is the loading matrix, then we seek B(find(fixed.position)) = fixed.value. Therefore, fixed.value must be a matrix of the same size as the loadings matrix and with the corresponding elements to be fixed at their appropriate values. All other elements of fixed.value are disregarded.
 
* '''fixed.weight''': [ ], a scalar (0 ? fixed.weight ? 1) indicating how strongly the fixed.value is imposed. A value of 0 (zero) does not impose the constraint at all, whereas a value of 1 (one) fixes the constraint.
 
* '''ridge.weight''': [ ], a scalar value between 0 and 1 that introduces a ridging in the update of the loading matrix. It is a penalty on the size of the esimated loadings. The closer to 1, the higher the ridge. Ridging is useful when a problem is difficult to fit.
 
* '''equality.G''': [ ], matrix with ''N'' columns, where ''N'' is the number of factors, used with equality.H. If '''A''' is the loadings for this mode then the constraint is imposed such that '''GA'''<sup>T</sup> = '''H'''. For example, if '''G''' is a row vector of ones and '''H''' is a vector of ones (1's), this would impose closure.
 
* '''equality.H''': [ ], matrix of size consistent with the constriant imposed by equality.G.
 
* '''equality.weight''': [ ], a scalar (0 ? equality.weight ? 1) indicating how strongly the equality.H and equality.G is imposed. A value of 0 (zero) does not impose the constraint at all, whereas a value of 1 (one) fixes the constraint.
 
* '''leftprod''': [0], If the loading matrix, '''B''' is of size JxR, the leftprod is a matrx '''G''' of size JxM. The loading '''B''' is then constrained to be of the form '''B = GH''', where only '''H''' is updated. For example, '''G''' may be a certain JxJ subspace, if the loadings are to be within a certain subspace.
 
* '''rightprod''': [0], If the loading matrix, '''B''' is of size JxR, the rightprod is a matrx '''G''' of size MxR. The loading '''B''' is then constrained to be of the form '''B''' = '''HG''', where only '''H''' is updated. For example, if rightprod is [1 1 0;0 0 1], then the first two components in '''B''' are forced to be the same.


* '''iterate_to_conv''': [0], Usually the constraints are imposed within an iterative algorithm. Some of the constraints use iterative algorithms themselves. Setting iterate_to_conv to one, will force the iterative constraint algorithms to continue until convergence.
To decompose a matrix (x) without non-negativity constraints use:


* '''timeaxis''': [], This field (if supplied) is used as the time axis when fitting loadings to a function (e.g. see exponential). Therefore, it must have the same number of elements as one of the loading vectors for this mode.
:options = als('options');


* '''description''': [1x1592 char],
:options.ccon = 'none';


If the constraint in a mode is set as fixed, then the loadings of that mode will not be updated, hence the initial loadings stay fixed.
:options.scon = 'none';
 
===Examples===


parafac demo gives a demonstration of the use of the PARAFAC algorithm.
:[c,s] = als(x,c0,options);


model = parafac(X,5) fits a five-component PARAFAC model to the array X using default settings.
The following shows an example of using soft-constraints on the second spectral component of a three-component solution assuming that the variable softs contains the spectrum to which component two should be constrained.


pred = parafac(Z,model) fits a parafac model to new data Z. The scores will be taken to be in the first mode, but you can change this by setting options.samplemodex to the mode which is the sample mode. Note, that the sample-mode dimension may be different for the old model and the new data, but all other dimensions must be the same.
:[m,n] = size(x);


options = parafac('options'); generates a set of default settings for PARAFAC. options.plots = 0; sets the plotting off.
:options = als('options');


options.init = 3; sets the initialization of PARAFAC to orthogonalized random numbers.
:options.sc = NaN\*ones(3,n);   %all 3 unconstrained


options.samplemodex = 2; Defines the second mode to be the sample-mode. Useful, for example, when fitting an existing model to new data has to provide the scores in the second mode.
:options.sc(2,:) = softs;      %constrain component 2


model = parafac(X,2,options); fits a two-component PARAFAC model with the settings defined in options.
:options.scwts = 0.5;          %consider as 1/2 of total signal in X


parafac io shows the I/O of the algorithm.
:[c,s] = als(x,c0,options);


===See Also===
===See Also===


[[datahat]], [[explode]], [[gram]], [[mpca]], [[outerm]], [[parafac2]], [[tld]], [[tucker]], [[unfoldm]]
[[mcr]], [[parafac]], [[pca]]

Revision as of 11:46, 25 September 2008

Purpose

Alternating Least Squares computational engine for multivariate curve resolution (MCR).

Synopsis

[c,s] = als(x,c0,options);

Description

ALS decomposes a matrix X as CS such that X = CS + E where E is minimized in a least squares sense.

Inputs are the matrix to be decomposed (x) (M by N), and the initial guess (c0). If (c0) is size M by K, where K is the number of factors, then it is assumed to be the initial guess for C. If (c0) is size K by N then it is assumed to be the initial guess for S. If M = K then, (c0) is assumed to be the initial guess for C.

An optional input (options) is described below.

The outputs are the estimated matrices C, (c), (M by K) and S, (s), (K by N). For a typical decomposition, (c) is a matrix of contributions and (s) is a matrix of spectra. The function

[c,s] = als(x,c0)

will decompose (x) using an non-negatively constrained alternating least squares. To include other constraints, use the (options) input described below.

Note that if no non-zero equality constraints are imposed on a factor the spectra are normalized to unit length. This can lead to significant scaling differences between factors that have non-zero equality constraints and those that do not.

Options

  • display: [ 'off' | {'on'} ] governs level of display to command window.
  • plots: [ 'none' | {'final'} ] governs level of plotting.
  • ccon: [ 'none' | 'reset' | {'fastnnls'} ] non-negativity on contributions. (fastnnls = true least-squares solution)
  • scon: [ 'none' | 'reset' | {'fastnnls'} ] non-negativity on spectra. (fastnnls = true least-squares solution)
  • cc: [ ] contributions equality constraints, must be a matrix with M rows and up to K columns with NaN where equality constraints are not applied and real value of the constraint where they are applied. If fewer than K columns are supplied, the missing columns will be filled in as unconstrained,
  • ccwts: [inf] a scalar value or a 1xK vector with elements corresponding to weightings on constraints (0, no constraint, 0<wt<inf imposes constraint "softly", and inf is hard constrained). If a scalar value is passed for (ccwts), that value is applied for all K factors.
Soft constraints are imposed by augmenting an additional column onto X and a weight of 1 weights a given equality constraint so that it has an equal influence as the average single variable in X. For additional information see P.J. Gemperline and E. Cash, “Advantages of Soft versus Hard Constraints in Self–Modeling Curve Resolution Problems. Alternating Least Squares with Penalty Functions”, Anal. Chem., 75(16), 4236–4243 (2003).
  • sc: [ ] spectra equality constraints, must be a matrix with N columns and up to K rows with NaN where equality contraints are not applied and real value of the constraint where they are applied. If fewer than K rows are supplied, the missing rows will be filled in as unconstrained.
  • scwts: [inf] weighting for spectral equality constraints (see ccwts) above.
  • sclc: [ ] contributions scale axis, vector with M elements otherwise 1:M is used.
  • scls: [ ] spectra scale axis, vector with N elements otherwise 1:N is used.
  • condition: [{'none'}| 'norm' ] type of conditioning to perform on S and C before each regression step. 'norm' uses conditioning which can help stabilize the regression when factors are significantly different in magnitude.
  • tolc: [ {1e-5} ] tolerance on non-negativity for contributions.
  • tols: [ {1e-5} ] tolerance on non-negativity for spectra.
  • ittol: [ {1e-8} ] convergence tolerance.
  • itmax: [ {100} ] maximum number of iterations.
  • timemax: [ {3600} ] maximum time for iterations.
  • rankfail: [ 'drop' |{'reset'}| 'random' | 'fail' ] how are rank deficiencies handled:
  • drop - drop deficient components from model,
  • reset - reset deficient components to initial guess,
  • random - replace deficient components with random vector,
  • fail - stop analysis, give error.

Examples

To decompose a matrix (x) without non-negativity constraints use:

options = als('options');
options.ccon = 'none';
options.scon = 'none';
[c,s] = als(x,c0,options);

The following shows an example of using soft-constraints on the second spectral component of a three-component solution assuming that the variable softs contains the spectrum to which component two should be constrained.

[m,n] = size(x);
options = als('options');
options.sc = NaN\*ones(3,n); %all 3 unconstrained
options.sc(2,:) = softs; %constrain component 2
options.scwts = 0.5; %consider as 1/2 of total signal in X
[c,s] = als(x,c0,options);

See Also

mcr, parafac, pca