Decomposition¶

fdatools.decomposition: functional principal component analysis (FPCA) and functional canonical correlation analysis (FCCA), both scikit-learn estimators.

Functional principal components and canonical correlation.

Two estimators replace R fda's pca.fd, varmx.pca.fd and cca.fd. Both work on the basis coefficients of an :class:~fdatools.core.FData and both follow the scikit-learn estimator API, so they drop straight into a :class:~sklearn.pipeline.Pipeline behind :class:~fdatools.smoothing.Smoother.

Functional PCA with a roughness penalty (Ramsay & Silverman §9.4) maximises the sample variance of the scores :math:\int \xi(t) x_i(t)\,dt while charging the harmonic :math:\xi for its own roughness. Writing :math:W = \int \varphi\varphi^{T} for the Gram matrix of the basis, :math:R = \int (L\varphi)(L\varphi)^{T} for the roughness penalty and :math:V = C C^{T}/N for the sample covariance of the centred coefficients, the harmonic coefficients solve the symmetric generalised eigenproblem

.. math::

W V W\,b = \mu\,(W + \lambda R)\,b ,
\qquad b^{T}(W + \lambda R)\,b = 1 .

With :math:\lambda = 0 this collapses to the unpenalised problem :math:V W b = \mu b, i.e. the eigen-decomposition of the covariance operator.

Functional CCA (chapter 11) is the same construction on the cross-covariance of two samples: with :math:S_{xy} = W_1 C_x C_y^{T} W_2 / N the canonical weights solve

.. math::

\begin{pmatrix} 0 & S_{xy} \\ S_{yx} & 0 \end{pmatrix} v
= \rho \begin{pmatrix} S_{xx} + \lambda_1 R_1 & 0 \\
                       0 & S_{yy} + \lambda_2 R_2 \end{pmatrix} v .

Every Gram matrix is exact (Gauss-Legendre on the basis breaks), where R's pca.fd uses a Romberg approximation good to four or five digits; see tests/parity/test_decomposition.py for the measured consequences.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline, FData
>>> from fdatools.decomposition import FPCA
>>> rng = np.random.default_rng(0)
>>> basis = BSpline(domain=(0.0, 1.0), n_basis=8)
>>> fd = FData(rng.standard_normal((8, 30)), basis)
>>> pca = FPCA(n=3).fit(fd)
>>> pca.scores.shape
(30, 3)
>>> bool(np.all(np.diff(pca.values) <= 1e-12))
True

FPCA ¶

FPCA(
    n: SupportsIndex = 2,
    *,
    lam: float | str = 0.0,
    penalty: int | LDO = 2,
    center: bool = True,
    basis: Basis | None = None,
)

Bases: TransformerMixin, BaseEstimator

Functional principal component analysis, optionally with a roughness penalty.

Replaces R's pca.fd and varmx.pca.fd.

Parameters:

Name	Type	Description	Default
`n`	`int or SupportsIndex`	Number of harmonics to keep. Default `2`. Any integer-like value (`numpy.int64`, ...) works; booleans and floats are rejected.	`2`
`lam`	`float or str`	Roughness penalty on the harmonics. `0.0` (default) is the unpenalised problem; `"gcv"` picks the penalty by leave-one-curve-out cross-validation of the `n`-harmonic reconstruction error over the log grid `10 -8, 10 -7, ..., 10 ** 8`.	`0.0`
`penalty`	`int or LDO`	Roughness operator `L`. Default `2` (the second derivative).	`2`
`center`	`bool`	Subtract the sample mean function before decomposing. Default `True`.	`True`
`basis`	`Basis`	Basis to attach to a plain coefficient matrix. Ignored when `fit` receives an :class:`~fdatools.core.FData`. `None` builds a cubic B-spline on `(0, 1)` with one basis function per column, which lets the estimator run inside a generic pipeline; pass the smoother's basis to decompose in the right metric.	`None`

Attributes:

Name	Type	Description
`harmonics`	`FData`	The `n` principal component functions.
`values`	`array`	The full eigenvalue spectrum, descending -- one entry per basis function (times the number of variables), as R reports it. After :meth:`rotate` it holds the variance of each rotated component instead.
`scores`	`array`	Score matrix of shape `(n_curves, n)`. For multivariate curves the per-variable contributions are summed; see :attr:`scores_by_var`.
`varprop`	`array`	Proportion of total variance explained by each retained harmonic.
`mean_fd`	`FData`	The sample mean function (a zero function when `center=False`).
`rotation`	`array or None`	The rotation applied by :meth:`rotate`, or `None` for an unrotated fit.

Notes

Eigenvector signs are mathematically arbitrary. fdatools makes the coefficient sum of every harmonic positive -- the rule R's pca.fd follows -- so repeated fits give identical output and match R's signs.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline, FData
>>> from fdatools.decomposition import FPCA
>>> rng = np.random.default_rng(1)
>>> fd = FData(rng.standard_normal((7, 40)), BSpline(domain=(0.0, 1.0), n_basis=7))
>>> pca = FPCA(n=2).fit(fd)
>>> pca.harmonics.n_curves
2
>>> bool(pca.varprop.sum() <= 1.0 + 1e-12)
True

harmonics `property` ¶

harmonics: FData

The principal component functions.

values `property` ¶

values: Array

The full eigenvalue spectrum, descending.

scores `property` ¶

scores: Array

Score matrix of shape (n_curves, n).

scores_by_var `property` ¶

scores_by_var: Array

Score contributions of shape (n_curves, n, n_vars).

varprop `property` ¶

varprop: Array

Proportion of total variance carried by each retained harmonic.

mean_fd `property` ¶

mean_fd: FData

The sample mean function.

rotation `property` ¶

rotation: Array | None

The rotation applied by :meth:rotate, or None.

fit ¶

fit(X: Any, y: Any = None) -> FPCA

Decompose X into functional principal components.

Parameters:

Name	Type	Description	Default
`X`	`FData or array of shape (n_samples, n_basis)`	Curves, or their basis coefficients one curve per row.	required
`y`	`ignored`	Present for API compatibility.	`None`

Returns:

Type	Description
`FPCA`	The fitted estimator.

Raises:

Type	Description
`ValueError`	If `n` is not a positive integer or `lam` is invalid.

Examples:

>>> import numpy as np
>>> from fdatools import Fourier, FData
>>> from fdatools.decomposition import FPCA
>>> rng = np.random.default_rng(3)
>>> fd = FData(rng.standard_normal((5, 20)), Fourier(domain=(0.0, 1.0), n_basis=5))
>>> FPCA(n=2).fit(fd).harmonics.n_curves
2

transform ¶

transform(X: Any) -> Array

Return the principal component scores of X.

Parameters:

Name	Type	Description	Default
`X`	`FData or array of shape (n_samples, n_basis)`	Curves to project onto the fitted harmonics.	required

Returns:

Type	Description
`array`	Scores of shape `(n_samples, n)`, summed over variables.

Raises:

Type	Description
`ValueError`	If `X` does not carry the number of variables seen during `fit`.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline, FData
>>> from fdatools.decomposition import FPCA
>>> rng = np.random.default_rng(4)
>>> fd = FData(rng.standard_normal((6, 25)), BSpline(domain=(0.0, 1.0), n_basis=6))
>>> pca = FPCA(n=2).fit(fd)
>>> bool(np.allclose(pca.transform(fd), pca.scores))
True

inverse_transform ¶

inverse_transform(X: Any) -> FData

Rebuild curves from scores.

Parameters:

Name	Type	Description	Default
`X`	`array of shape (n_samples, n)`	Principal component scores.	required

Returns:

Type	Description
`FData`	`mean_fd` plus the score-weighted harmonics.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline, FData
>>> from fdatools.decomposition import FPCA
>>> rng = np.random.default_rng(5)
>>> fd = FData(rng.standard_normal((4, 30)), BSpline(domain=(0.0, 1.0), n_basis=4))
>>> pca = FPCA(n=4).fit(fd)
>>> bool(np.allclose(pca.inverse_transform(pca.scores).coefs, fd.coefs))
True

rotate ¶

rotate(method: str = 'varimax') -> FPCA

Return a copy of this fit with rotated harmonics.

Parameters:

Name	Type	Description	Default
`method`	`str`	Only `"varimax"` is available; it maximises the variance of the squared harmonic values on a 501-point grid, exactly the quantity R's `varmx.pca.fd` rotates.	`'varimax'`

Returns:

Type	Description
`FPCA`	A new fitted estimator whose `harmonics`, `scores`, `values` and `varprop` describe the rotated solution: `values` is the variance of each rotated component and `varprop` its share of the total variance of the unrotated spectrum. `mean_fd` is unchanged. Rotated harmonics keep the positive-coefficient-sum sign rule.

Raises:

Type	Description
`ValueError`	If `method` is not `"varimax"`.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline, FData
>>> from fdatools.decomposition import FPCA
>>> rng = np.random.default_rng(6)
>>> fd = FData(rng.standard_normal((8, 30)), BSpline(domain=(0.0, 1.0), n_basis=8))
>>> rotated = FPCA(n=3).fit(fd).rotate("varimax")
>>> bool(np.allclose(rotated.rotation.T @ rotated.rotation, np.eye(3)))
True

plot ¶

plot(ax: Any = None, *, n_points: int = 201) -> list[Any]

Plot each harmonic as a perturbation of the mean.

Every panel shows the mean function and the mean plus and minus the harmonic scaled by the standard deviation of its scores -- the display R's plot.pca.fd produces.

Parameters:

Name	Type	Description	Default
`ax`	`matplotlib.axes.Axes or sequence of Axes`	Axes to draw on, one per harmonic. `None` creates one row of panels.	`None`
`n_points`	`int`	Number of evaluation points. Default `201`.	`201`

Returns:

Type	Description
`list of matplotlib.axes.Axes`	The axes drawn on.

Raises:

Type	Description
`ValueError`	If fewer axes than harmonics are supplied.

Examples:

>>> import matplotlib
>>> matplotlib.use("Agg")
>>> import numpy as np
>>> from fdatools import BSpline, FData
>>> from fdatools.decomposition import FPCA
>>> rng = np.random.default_rng(7)
>>> fd = FData(rng.standard_normal((6, 20)), BSpline(domain=(0.0, 1.0), n_basis=6))
>>> len(FPCA(n=2).fit(fd).plot())
2

FCCA ¶

FCCA(
    n: SupportsIndex = 2,
    *,
    lam1: float = 0.0,
    lam2: float = 0.0,
    penalty: int | LDO = 2,
    center: bool = True,
)

Bases: BaseEstimator

Functional canonical correlation analysis.

Replaces R's cca.fd.

Parameters:

Name	Type	Description	Default
`n`	`int or SupportsIndex`	Number of canonical variate pairs to keep. Default `2`. Any integer-like value (`numpy.int64`, ...) works; booleans and floats are rejected.	`2`
`lam1`	`float`	Roughness penalties on the first and second set of canonical weights. Without a penalty the problem is degenerate whenever there are fewer curves than basis functions (every correlation is one).	`0.0`
`lam2`	`float`	Roughness penalties on the first and second set of canonical weights. Without a penalty the problem is degenerate whenever there are fewer curves than basis functions (every correlation is one).	`0.0`
`penalty`	`int or LDO`	Roughness operator `L`. Default `2`.	`2`
`center`	`bool`	Subtract each sample mean before decomposing. Default `True`.	`True`

Attributes:

Name	Type	Description
`weights1, weights2`	`FData`	Canonical weight functions, normalised to unit :math:`L^2` norm. The sign of each pair is fixed so the first weight function has a positive coefficient sum; both members of a pair share that sign, so their scores stay positively correlated.
`correlations`	`array`	The full canonical correlation spectrum, descending, as R reports it.
`scores1, scores2`	`array`	Canonical variate scores of shape `(n_curves, n)`.

Examples:

>>> import numpy as np
>>> from fdatools import BSpline, FData
>>> from fdatools.decomposition import FCCA
>>> rng = np.random.default_rng(2)
>>> basis = BSpline(domain=(0.0, 1.0), n_basis=6)
>>> a = rng.standard_normal((6, 25))
>>> fd1 = FData(a, basis)
>>> fd2 = FData(a + 0.1 * rng.standard_normal((6, 25)), basis)
>>> cca = FCCA(n=2, lam1=1e-6, lam2=1e-6).fit(fd1, fd2)
>>> bool(cca.correlations[0] > 0.9)
True

weights1 `property` ¶

weights1: FData

Canonical weight functions of the first sample.

weights2 `property` ¶

weights2: FData

Canonical weight functions of the second sample.

correlations `property` ¶

correlations: Array

The full canonical correlation spectrum, descending.

scores1 `property` ¶

scores1: Array

Canonical variate scores of the first sample.

scores2 `property` ¶

scores2: Array

Canonical variate scores of the second sample.

fit ¶

fit(X: FData, y: FData) -> FCCA

Find the canonical weight functions of two samples of curves.

Parameters:

Name	Type	Description	Default
`X`	`FData`	Two samples observed on the same curves, in that order.	required
`y`	`FData`	Two samples observed on the same curves, in that order.	required

Returns:

Type	Description
`FCCA`	The fitted estimator.

Raises:

Type	Description
`TypeError`	If either argument is not an :class:`~fdatools.core.FData`.
`ValueError`	If the samples hold different numbers of curves, are multivariate, or a parameter is out of range.

Examples:

>>> import numpy as np
>>> from fdatools import Fourier, FData
>>> from fdatools.decomposition import FCCA
>>> rng = np.random.default_rng(8)
>>> basis = Fourier(domain=(0.0, 1.0), n_basis=5)
>>> x = rng.standard_normal((5, 40))
>>> cca = FCCA(n=1, lam1=1e-3, lam2=1e-3).fit(FData(x, basis), FData(-x, basis))
>>> cca.scores1.shape
(40, 1)

transform ¶

transform(X: FData, y: FData) -> tuple[Array, Array]

Return the canonical variate scores of two new samples.

Parameters:

Name	Type	Description	Default
`X`	`FData`	Curves to project, in the order given to :meth:`fit`.	required
`y`	`FData`	Curves to project, in the order given to :meth:`fit`.	required

Returns:

Type	Description
`tuple of array`	Scores of shape `(n_curves, n)` for each sample.

Raises:

Type	Description
`ValueError`	If a sample's basis size differs from the one seen during `fit`.

Examples:

>>> import numpy as np
>>> from fdatools import Fourier, FData
>>> from fdatools.decomposition import FCCA
>>> rng = np.random.default_rng(9)
>>> basis = Fourier(domain=(0.0, 1.0), n_basis=5)
>>> x = FData(rng.standard_normal((5, 30)), basis)
>>> y = FData(rng.standard_normal((5, 30)), basis)
>>> cca = FCCA(n=2, lam1=1e-2, lam2=1e-2).fit(x, y)
>>> bool(np.allclose(cca.transform(x, y)[0], cca.scores1))
True

Decomposition¶

FPCA ¶

harmonics property ¶

values property ¶

scores property ¶

scores_by_var property ¶

varprop property ¶

mean_fd property ¶

rotation property ¶

fit ¶

transform ¶

inverse_transform ¶

rotate ¶

plot ¶

FCCA ¶

weights1 property ¶

weights2 property ¶

correlations property ¶

scores1 property ¶

scores2 property ¶

fit ¶

transform ¶

harmonics `property` ¶

values `property` ¶

scores `property` ¶

scores_by_var `property` ¶

varprop `property` ¶

mean_fd `property` ¶

rotation `property` ¶

weights1 `property` ¶

weights2 `property` ¶

correlations `property` ¶

scores1 `property` ¶

scores2 `property` ¶