Hybrid Mapping Techniques

PCovR

class skmatter.decomposition.PCovR(mixing=0.5, n_components=None, svd_solver='auto', tol=1e-12, space='auto', regressor=None, iterated_power='auto', random_state=None, whiten=False)

Bases: RegressorMixin, MultiOutputMixin, _BasePCov

Principal Covariates Regression (PCovR).

As described in [deJong1992], PCovR determines a latent-space projection \(\mathbf{T}\) which minimizes a combined loss in supervised and unsupervised tasks.

This projection is determined by the eigendecomposition of a modified gram matrix \(\mathbf{\tilde{K}}\)

\[\mathbf{\tilde{K}} = \alpha \mathbf{X} \mathbf{X}^T + (1 - \alpha) \mathbf{\hat{Y}}\mathbf{\hat{Y}}^T\]

where \(\alpha\) is a mixing parameter and \(\mathbf{X}\) and \(\mathbf{\hat{Y}}\) are matrices of shapes \((n_{samples}, n_{features})\) and \((n_{samples}, n_{properties})\), respectively, which contain the input and approximate targets. For \((n_{samples} < n_{features})\), this can be more efficiently computed using the eigendecomposition of a modified covariance matrix \(\mathbf{\tilde{C}}\)

\[\mathbf{\tilde{C}} = \alpha \mathbf{X}^T \mathbf{X} + (1 - \alpha) \left(\left(\mathbf{X}^T \mathbf{X}\right)^{-\frac{1}{2}} \mathbf{X}^T \mathbf{\hat{Y}}\mathbf{\hat{Y}}^T \mathbf{X} \left(\mathbf{X}^T \mathbf{X}\right)^{-\frac{1}{2}}\right)\]

For all PCovR methods, it is strongly suggested that \(\mathbf{X}\) and \(\mathbf{Y}\) are centered and scaled to unit variance, otherwise the results will change drastically near \(\alpha \to 0\) and \(\alpha \to 1\). This can be done with the companion preprocessing classes, where

>>> from skmatter.preprocessing import StandardFlexibleScaler as SFS
>>> import numpy as np
>>>
>>> # Set column_wise to True when the columns are relative to one another,
>>> # False otherwise.
>>> scaler = SFS(column_wise=True)
>>>
>>> A = np.array([[1, 2], [2, 1]])  # replace with your matrix
>>> scaler.fit(A)
StandardFlexibleScaler(column_wise=True)
>>> A = scaler.transform(A)

Parameters:

• mixing (float, default=0.5) – mixing parameter, as described in PCovR as \({\alpha}\), here named to avoid confusion with regularization parameter alpha

• n_components (int, float or str, default=None) –

  Number of components to keep. if n_components is not set all components are kept:

  n_components == min(n_samples, n_features)
  

• svd_solver ({'auto', 'full', 'arpack', 'randomized'}, default='auto') –

  If auto :

  The solver is selected by a default policy based on X.shape and n_components: if the input data is larger than 500x500 and the number of components to extract is lower than 80% of the smallest dimension of the data, then the more efficient ‘randomized’ method is enabled. Otherwise the exact full SVD is computed and optionally truncated afterwards.

  If full :

  run exact full SVD calling the standard LAPACK solver via scipy.linalg.svd and select the components by postprocessing

  If arpack :

  run SVD truncated to n_components calling ARPACK solver via scipy.sparse.linalg.svds. It requires strictly 0 < n_components < min(X.shape)

  If randomized :

  run randomized SVD by the method of Halko et al.

• tol (float, default=1e-12) – Tolerance for singular values computed by svd_solver == ‘arpack’. Must be of range [0.0, infinity).

• space ({'feature', 'sample', 'auto'}, default='auto') – whether to compute the PCovC in sample or feature space. The default is equal to sample when \({n_{samples} < n_{features}}\) and feature when \({n_{features} < n_{samples}}\)

• regressor ({Ridge, RidgeCV, LinearRegression, precomputed}, default=None) – regressor for computing approximated \({\mathbf{\hat{Y}}}\). The regressor should be one sklearn.linear_model.Ridge, sklearn.linear_model.RidgeCV, or sklearn.linear_model.LinearRegression. If a pre-fitted regressor is provided, it is used to compute \({\mathbf{\hat{Y}}}\). Note that any pre-fitting of the regressor will be lost if PCovR is within a composite estimator that enforces cloning, e.g., sklearn.compose.TransformedTargetRegressor or sklearn.pipeline.Pipeline with model caching. In such cases, the regressor will be re-fitted on the same training data as the composite estimator. If precomputed, we assume that the y passed to the fit function is the regressed form of the targets \({\mathbf{\hat{Y}}}\). If None, sklearn.linear_model.Ridge('alpha':1e-6, 'fit_intercept':False, 'tol':1e-12) is used as the regressor.

• iterated_power (int or 'auto', default='auto') – Number of iterations for the power method computed by svd_solver == ‘randomized’. Must be of range [0, infinity).

• random_state (int, numpy.random.RandomState instance or None, default=None) – Used when the ‘arpack’ or ‘randomized’ solvers are used. Pass an int for reproducible results across multiple function calls.

• whiten (boolean, deprecated)

mixing

mixing parameter, as described in PCovR as \({\alpha}\)

Type:

float, default=0.5

tol

Tolerance for singular values computed by svd_solver == ‘arpack’. Must be of range [0.0, infinity).

Type:

float, default=1e-12

space

whether to compute the PCovR in sample or feature space. The default is equal to sample when \({n_{samples} < n_{features}}\) and feature when \({n_{features} < n_{samples}}\)

Type:

{‘feature’, ‘sample’, ‘auto’}, default=’auto’

n_components_

The estimated number of components, which equals the parameter n_components, or the lesser value of n_features and n_samples if n_components is None.

Type:

int

pxt_

the projector, or weights, from the input space \(\mathbf{X}\) to the latent-space projection \(\mathbf{T}\)

Type:

numpy.ndarray of size \(({n_{samples}, n_{components}})\)

pxy_

the projector, or weights, from the input space \(\mathbf{X}\) to the properties \(\mathbf{Y}\)

Type:

numpy.ndarray of size \(({n_{samples}, n_{properties}})\)

pty_

the projector, or weights, from the latent-space projection \(\mathbf{T}\) to the properties \(\mathbf{Y}\)

Type:

numpy.ndarray of size \(({n_{components}, n_{properties}})\)

explained_variance_

The amount of variance explained by each of the selected components. Equal to n_components largest eigenvalues of the PCovR-modified covariance matrix of \(\mathbf{X}\).

Type:

numpy.ndarray of shape (n_components,)

singular_values_

The singular values corresponding to each of the selected components.

Type:

numpy.ndarray of shape (n_components,)

Examples

>>> import numpy as np
>>> from skmatter.decomposition import PCovR
>>> X = np.array([[-1, 1, -3, 1], [1, -2, 1, 2], [-2, 0, -2, -2], [1, 0, 2, -1]])
>>> Y = np.array([[0, -5], [-1, 1], [1, -5], [-3, 2]])
>>> pcovr = PCovR(mixing=0.1, n_components=2)
>>> pcovr.fit(X, Y)
PCovR(mixing=0.1, n_components=2)
>>> pcovr.transform(X)
array([[ 3.2630561 ,  0.06663787],
       [-2.69395511, -0.41582771],
       [ 3.48683147, -0.83164387],
       [-4.05593245,  1.18083371]])
>>> pcovr.predict(X)
array([[ 0.01371776, -5.00945512],
       [-1.02805338,  1.06736871],
       [ 0.98166504, -4.98307078],
       [-2.9963189 ,  1.98238856]])

fit(X, Y, W=None)

Fit the model with X and Y. Depending on the dimensions of X, calls either _fit_feature_space or _fit_sample_space

Parameters:

• X (numpy.ndarray, shape (n_samples, n_features)) –

  Training data, where n_samples is the number of samples and n_features is the number of features.

  It is suggested that \(\mathbf{X}\) be centered by its column- means and scaled. If features are related, the matrix should be scaled to have unit variance, otherwise \(\mathbf{X}\) should be scaled so that each feature has a variance of 1 / n_features.

• Y (numpy.ndarray, shape (n_samples, n_properties)) –

  Training data, where n_samples is the number of samples and n_properties is the number of properties

  It is suggested that \(\mathbf{Y}\) be centered by its column- means and scaled. If features are related, the matrix should be scaled to have unit variance, otherwise \(\mathbf{Y}\) should be scaled so that each feature has a variance of 1 / n_features.

  If the passed regressor = precomputed, it is assumed that Y is the regressed form of the properties, \({\mathbf{\hat{Y}}}\).

• W (numpy.ndarray, shape (n_features, n_properties)) – Regression weights, optional when regressor is precomputed. If not passed, it is assumed that W = np.linalg.lstsq(X, Y, self.tol)[0]

_fit_feature_space(X, Y, Yhat)

In feature-space PCovR, the projectors are determined by:

\[\mathbf{\tilde{C}} = \alpha \mathbf{X}^T \mathbf{X} + (1 - \alpha) \left(\left(\mathbf{X}^T \mathbf{X}\right)^{-\frac{1}{2}} \mathbf{X}^T \mathbf{\hat{Y}}\mathbf{\hat{Y}}^T \mathbf{X} \left(\mathbf{X}^T \mathbf{X}\right)^{-\frac{1}{2}}\right)\]

where

\[\mathbf{P}_{XT} = (\mathbf{X}^T \mathbf{X})^{-\frac{1}{2}} \mathbf{U}_\mathbf{\tilde{C}}^T \mathbf{\Lambda}_\mathbf{\tilde{C}}^{\frac{1}{2}}\]

\[\mathbf{P}_{TX} = \mathbf{\Lambda}_\mathbf{\tilde{C}}^{-\frac{1}{2}} \mathbf{U}_\mathbf{\tilde{C}}^T (\mathbf{X}^T \mathbf{X})^{\frac{1}{2}}\]

\[\mathbf{P}_{TY} = \mathbf{\Lambda}_\mathbf{\tilde{C}}^{-\frac{1}{2}} \mathbf{U}_\mathbf{\tilde{C}}^T (\mathbf{X}^T \mathbf{X})^{-\frac{1}{2}} \mathbf{X}^T \mathbf{Y}\]

_fit_sample_space(X, Y, Yhat, W)

In sample-space PCovR, the projectors are determined by:

\[\mathbf{\tilde{K}} = \alpha \mathbf{X} \mathbf{X}^T + (1 - \alpha) \mathbf{\hat{Y}}\mathbf{\hat{Y}}^T\]

where

\[\mathbf{P}_{XT} = \left(\alpha \mathbf{X}^T + (1 - \alpha) \mathbf{W} \mathbf{\hat{Y}}^T\right) \mathbf{U}_\mathbf{\tilde{K}} \mathbf{\Lambda}_\mathbf{\tilde{K}}^{-\frac{1}{2}}\]

\[\mathbf{P}_{TX} = \mathbf{\Lambda}_\mathbf{\tilde{K}}^{-\frac{1}{2}} \mathbf{U}_\mathbf{\tilde{K}}^T \mathbf{X}\]

\[\mathbf{P}_{TY} = \mathbf{\Lambda}_\mathbf{\tilde{K}}^{-\frac{1}{2}} \mathbf{U}_\mathbf{\tilde{K}}^T \mathbf{Y}\]

transform(X=None)

Apply dimensionality reduction to X.

X is projected on the first principal components as determined by the modified PCovR distances.

Parameters:

X (numpy.ndarray, shape (n_samples, n_features)) – New data, where n_samples is the number of samples and n_features is the number of features.

predict(X=None, T=None)

Predicts the property values using regression on X or T.

inverse_transform(T)

Transform data back to its original space.

\[\mathbf{\hat{X}} = \mathbf{T} \mathbf{P}_{TX} = \mathbf{X} \mathbf{P}_{XT} \mathbf{P}_{TX}\]

Parameters:

T (ndarray, shape (n_samples, n_components)) – Projected data, where n_samples is the number of samples and n_components is the number of components.

Returns:

X_original (numpy.ndarray, shape (n_samples, n_features))

score(X, y, T=None)

Return the (negative) total reconstruction error for X and Y, defined as:

\[\ell_{X} = \frac{\lVert \mathbf{X} - \mathbf{T}\mathbf{P}_{TX} \rVert ^ 2} {\lVert \mathbf{X}\rVert ^ 2}\]

and

\[\ell_{Y} = \frac{\lVert \mathbf{Y} - \mathbf{T}\mathbf{P}_{TY} \rVert ^ 2} {\lVert \mathbf{Y}\rVert ^ 2}\]

The negative loss \(-\ell = -(\ell_{X} + \ell{Y})\) is returned for easier use in sklearn pipelines, e.g., a grid search, where methods named ‘score’ are meant to be maximized.

Parameters:

• X (numpy.ndarray of shape (n_samples, n_features)) – The data.

• Y (numpy.ndarray of shape (n_samples, n_properties)) – The target.

Returns:

loss (float) – Negative sum of the loss in reconstructing X from the latent-space projection T and the loss in predicting Y from the latent-space projection T

PCovC

class skmatter.decomposition.PCovC(mixing=0.5, n_components=None, svd_solver='auto', tol=1e-12, z_mean_tol=1e-12, z_var_tol=1.5, space='auto', classifier=None, scale_z=False, iterated_power='auto', random_state=None, whiten=False)

Bases: LinearClassifierMixin, _BasePCov

Principal Covariates Classification (PCovC).

As described in [Jorgensen2025], PCovC determines a latent-space projection \(\mathbf{T}\) which minimizes a combined loss in supervised and unsupervised tasks.

This projection is determined by the eigendecomposition of a modified gram matrix \(\mathbf{\tilde{K}}\)

\[\mathbf{\tilde{K}} = \alpha \mathbf{X} \mathbf{X}^T + (1 - \alpha) \mathbf{Z}\mathbf{Z}^T\]

where \(\alpha\) is a mixing parameter, \(\mathbf{X}\) is an input matrix of shape \((n_{samples}, n_{features})\), and \(\mathbf{Z}\) is a matrix of class confidence scores of shape \((n_{samples}, n_{classes})\). For \((n_{samples} < n_{features})\), this can be more efficiently computed using the eigendecomposition of a modified covariance matrix \(\mathbf{\tilde{C}}\)

\[\mathbf{\tilde{C}} = \alpha \mathbf{X}^T \mathbf{X} + (1 - \alpha) \left(\left(\mathbf{X}^T \mathbf{X}\right)^{-\frac{1}{2}} \mathbf{X}^T \mathbf{Z}\mathbf{Z}^T \mathbf{X} \left(\mathbf{X}^T \mathbf{X}\right)^{-\frac{1}{2}}\right)\]

For all PCovC methods, it is strongly suggested that \(\mathbf{X}\) is centered and scaled to unit variance, otherwise the results will change drastically near \(\alpha \to 0\) and \(\alpha \to 1\). This can be done with the companion preprocessing classes, where

>>> from skmatter.preprocessing import StandardFlexibleScaler as SFS
>>> import numpy as np
>>>
>>> # Set column_wise to True when the columns are relative to one another,
>>> # False otherwise.
>>> scaler = SFS(column_wise=True)
>>>
>>> A = np.array([[1, 2], [2, 1]])  # replace with your matrix
>>> scaler.fit(A)
StandardFlexibleScaler(column_wise=True)
>>> A = scaler.transform(A)

Parameters:

• mixing (float, default=0.5) – mixing parameter, as described in PCovC as \({\alpha}\), here named to avoid confusion with regularization parameter alpha

• n_components (int, float or str, default=None) –

  Number of components to keep. if n_components is not set all components are kept:

  n_components == min(n_samples, n_features)
  

• svd_solver ({'auto', 'full', 'arpack', 'randomized'}, default='auto') –

  If auto :

  The solver is selected by a default policy based on X.shape and n_components: if the input data is larger than 500x500 and the number of components to extract is lower than 80% of the smallest dimension of the data, then the more efficient ‘randomized’ method is enabled. Otherwise the exact full SVD is computed and optionally truncated afterwards.

  If full :

  run exact full SVD calling the standard LAPACK solver via scipy.linalg.svd and select the components by postprocessing

  If arpack :

  run SVD truncated to n_components calling ARPACK solver via scipy.sparse.linalg.svds. It requires strictly 0 < n_components < min(X.shape)

  If randomized :

  run randomized SVD by the method of Halko et al.

• tol (float, default=1e-12) – Tolerance for singular values computed by svd_solver == ‘arpack’. Must be of range [0.0, infinity).

• z_mean_tol (float, default=1e-12) – Tolerance for the column means of Z. Must be of range [0.0, infinity).

• z_var_tol (float, default=1.5) – Tolerance for the column variances of Z. Must be of range [0.0, infinity).

• space ({'feature', 'sample', 'auto'}, default='auto') – whether to compute the PCovC in sample or feature space. The default is equal to sample when \({n_{samples} < n_{features}}\) and feature when \({n_{features} < n_{samples}}\)

• classifier (estimator object or precomputed, default=None) –

  classifier for computing \({\mathbf{Z}}\). The classifier should be one of the following:

  • sklearn.linear_model.LogisticRegression()

  • sklearn.linear_model.LogisticRegressionCV()

  • sklearn.svm.LinearSVC()

  • sklearn.discriminant_analysis.LinearDiscriminantAnalysis()

  • sklearn.linear_model.RidgeClassifier()

  • sklearn.linear_model.RidgeClassifierCV()

  • sklearn.linear_model.Perceptron()

  If a pre-fitted classifier is provided, it is used to compute \({\mathbf{Z}}\). Note that any pre-fitting of the classifier will be lost if PCovC is within a composite estimator that enforces cloning, e.g., sklearn.pipeline.Pipeline with model caching. In such cases, the classifier will be re-fitted on the same training data as the composite estimator. If None, sklearn.linear_model.LogisticRegression() is used as the classifier.

• scale_z (bool, default=False) – Whether to scale Z prior to eigendecomposition.

• iterated_power (int or 'auto', default='auto') – Number of iterations for the power method computed by svd_solver == ‘randomized’. Must be of range [0, infinity).

• random_state (int, RandomState instance or None, default=None) – Used when the ‘arpack’ or ‘randomized’ solvers are used. Pass an int for reproducible results across multiple function calls.

• whiten (boolean, deprecated)

mixing

mixing parameter, as described in PCovC as \({\alpha}\)

Type:

float, default=0.5

tol

Tolerance for singular values computed by svd_solver == ‘arpack’. Must be of range [0.0, infinity).

Type:

float, default=1e-12

z_mean_tol

Tolerance for the column means of Z. Must be of range [0.0, infinity).

Type:

float

z_var_tol

Tolerance for the column variances of Z. Must be of range [0.0, infinity).

Type:

float

space

whether to compute the PCovC in sample or feature space. The default is equal to sample when \({n_{samples} < n_{features}}\) and feature when \({n_{features} < n_{samples}}\)

Type:

{‘feature’, ‘sample’, ‘auto’}, default=’auto’

n_components_

The estimated number of components, which equals the parameter n_components, or the lesser value of n_features and n_samples if n_components is None.

Type:

int

classifier

The linear classifier passed for fitting.

Type:

estimator object

z_classifier_

The linear classifier fit between \(\mathbf{X}\) and \(\mathbf{Y}\).

Type:

estimator object

classifier_

The linear classifier fit between \(\mathbf{T}\) and \(\mathbf{Y}\).

Type:

estimator object

pxt_

the projector, or weights, from the input space \(\mathbf{X}\) to the latent-space projection \(\mathbf{T}\)

Type:

ndarray of size \(({n_{features}, n_{components}})\)

pxz_

the projector, or weights, from the input space \(\mathbf{X}\) to the class confidence scores \(\mathbf{Z}\)

Type:

ndarray of size \(({n_{features}, })\) or \(({n_{features}, n_{classes}})\)

ptz_

the projector, or weights, from the latent-space projection \(\mathbf{T}\) to the class confidence scores \(\mathbf{Z}\)

Type:

ndarray of size \(({n_{components}, })\) or \(({n_{components}, n_{classes}})\)

scale_z

Whether Z is being scaled prior to eigendecomposition

Type:

bool

explained_variance_

The amount of variance explained by each of the selected components. Equal to n_components largest eigenvalues of the PCovC-modified covariance matrix of \(\mathbf{X}\).

Type:

numpy.ndarray of shape (n_components,)

singular_values_

The singular values corresponding to each of the selected components.

Type:

numpy.ndarray of shape (n_components,)

Examples

>>> import numpy as np
>>> from skmatter.decomposition import PCovC
>>> from sklearn.preprocessing import StandardScaler
>>> X = np.array([[-1, 0, -2, 3], [3, -2, 0, 1], [-3, 0, -1, -1], [1, 3, 0, -2]])
>>> X = StandardScaler().fit_transform(X)
>>> Y = np.array([0, 1, 2, 0])
>>> pcovc = PCovC(mixing=0.1, n_components=2)
>>> pcovc.fit(X, Y)
PCovC(mixing=0.1, n_components=2)
>>> pcovc.transform(X)
array([[-0.4794854 , -0.46228114],
       [ 1.9416966 ,  0.2532831 ],
       [-1.08744947,  0.89117784],
       [-0.37476173, -0.6821798 ]])
>>> pcovc.predict(X)
array([0, 1, 2, 0])

fit(X, Y, W=None)

Fit the model with X and Y.

Note that W is taken from the coefficients of a linear classifier fit between X and Y to compute Z:

\[\mathbf{Z} = \mathbf{X} \mathbf{W}\]

We then call either _fit_feature_space or _fit_sample_space, using Z as our approximation of Y. Finally, we refit a classifier on T and Y to obtain \(\mathbf{P}_{TZ}\).

Parameters:

• X (numpy.ndarray, shape (n_samples, n_features)) –

  Training data, where n_samples is the number of samples and n_features is the number of features.

  It is suggested that \(\mathbf{X}\) be centered by its column- means and scaled. If features are related, the matrix should be scaled to have unit variance, otherwise \(\mathbf{X}\) should be scaled so that each feature has a variance of 1 / n_features.

• Y (numpy.ndarray, shape (n_samples,)) – Training data, where n_samples is the number of samples.

• W (numpy.ndarray, shape (n_features, n_classes)) – Classification weights, optional when classifier is precomputed. If not passed, it is assumed that the weights will be taken from a linear classifier fit between \(\mathbf{X}\) and \(\mathbf{Y}\)

_fit_feature_space(X, Y, Z)

In feature-space PCovC, the projectors are determined by:

\[\mathbf{\tilde{C}} = \alpha \mathbf{X}^T \mathbf{X} + (1 - \alpha) \left(\left(\mathbf{X}^T \mathbf{X}\right)^{-\frac{1}{2}} \mathbf{X}^T \mathbf{Z}\mathbf{Z}^T \mathbf{X} \left(\mathbf{X}^T \mathbf{X}\right)^{-\frac{1}{2}}\right)\]

where

\[\mathbf{P}_{XT} = (\mathbf{X}^T \mathbf{X})^{-\frac{1}{2}} \mathbf{U}_\mathbf{\tilde{C}}^T \mathbf{\Lambda}_\mathbf{\tilde{C}}^{\frac{1}{2}}\]

\[\mathbf{P}_{TX} = \mathbf{\Lambda}_\mathbf{\tilde{C}}^{-\frac{1}{2}} \mathbf{U}_\mathbf{\tilde{C}}^T (\mathbf{X}^T \mathbf{X})^{\frac{1}{2}}\]

_fit_sample_space(X, Y, Z, W)

In sample-space PCovC, the projectors are determined by:

\[\mathbf{\tilde{K}} = \alpha \mathbf{X} \mathbf{X}^T + (1 - \alpha) \mathbf{Z}\mathbf{Z}^T\]

where

\[\mathbf{P}_{XT} = \left(\alpha \mathbf{X}^T + (1 - \alpha) \mathbf{W} \mathbf{Z}^T\right) \mathbf{U}_\mathbf{\tilde{K}} \mathbf{\Lambda}_\mathbf{\tilde{K}}^{-\frac{1}{2}}\]

\[\mathbf{P}_{TX} = \mathbf{\Lambda}_\mathbf{\tilde{K}}^{-\frac{1}{2}} \mathbf{U}_\mathbf{\tilde{K}}^T \mathbf{X}\]

transform(X=None)

Apply dimensionality reduction to X.

X is projected on the first principal components as determined by the modified PCovC distances.

Parameters:

X (numpy.ndarray, shape (n_samples, n_features)) – New data, where n_samples is the number of samples and n_features is the number of features.

predict(X=None, T=None)

Predicts the property labels using classification on T.

inverse_transform(T)

Transform data back to its original space.

\[\mathbf{\hat{X}} = \mathbf{T} \mathbf{P}_{TX} = \mathbf{X} \mathbf{P}_{XT} \mathbf{P}_{TX}\]

Parameters:

T (ndarray, shape (n_samples, n_components)) – Projected data, where n_samples is the number of samples and n_components is the number of components.

Returns:

X_original (numpy.ndarray, shape (n_samples, n_features))

decision_function(X=None, T=None)

Predicts confidence scores from X or T.

\[\mathbf{Z} = \mathbf{T} \mathbf{P}_{TZ} = \mathbf{X} \mathbf{P}_{XT} \mathbf{P}_{TZ} = \mathbf{X} \mathbf{P}_{XZ}\]

Parameters:

• X (ndarray, shape(n_samples, n_features)) – Original data for which we want to get confidence scores, where n_samples is the number of samples and n_features is the number of features.

• T (ndarray, shape (n_samples, n_components)) – Projected data for which we want to get confidence scores, where n_samples is the number of samples and n_components is the number of components.

Returns:

Z (numpy.ndarray, shape (n_samples,) or (n_samples, n_classes)) – Confidence scores. For binary classification, has shape (n_samples,), for multiclass classification, has shape (n_samples, n_classes)

score(X, y, sample_weight=None)

Return accuracy on provided data and labels.

In multi-label classification, this is the subset accuracy which is a harsh metric since you require for each sample that each label set be correctly predicted.

Parameters:

• X (array-like of shape (n_samples, n_features)) – Test samples.

• y (array-like of shape (n_samples,) or (n_samples, n_outputs)) – True labels for X.

• sample_weight (array-like of shape (n_samples,), default=None) – Sample weights.

Returns:

score (float) – Mean accuracy of self.predict(X) w.r.t. y.

Kernel PCovR

class skmatter.decomposition.KernelPCovR(mixing=0.5, n_components=None, svd_solver='auto', regressor=None, kernel='linear', gamma=None, degree=3, coef0=1, kernel_params=None, center=False, fit_inverse_transform=False, tol=1e-12, n_jobs=None, iterated_power='auto', random_state=None)

Bases: _BaseKPCov

Kernel Principal Covariates Regression (KPCovR).

As described in [Helfrecht2020], KPCovR determines a latent-space projection \(\mathbf{T}\) which minimizes a combined loss in supervised and unsupervised tasks in the reproducing kernel Hilbert space (RKHS).

This projection is determined by the eigendecomposition of a modified gram matrix \(\mathbf{\tilde{K}}\)

\[\mathbf{\tilde{K}} = \alpha \mathbf{K} + (1 - \alpha) \mathbf{\hat{Y}}\mathbf{\hat{Y}}^T\]

where \(\alpha\) is a mixing parameter, \(\mathbf{K}\) is the input kernel of shape \((n_{samples}, n_{samples})\) and \(\mathbf{\hat{Y}}\) is the target matrix of shape \((n_{samples}, n_{properties})\).

Parameters:

• mixing (float, default=0.5) – mixing parameter, as described in PCovR as \({\alpha}\)

• n_components (int, float or str, default=None) –

  Number of components to keep. if n_components is not set all components are kept:

  n_components == n_samples
  

• svd_solver ({'auto', 'full', 'arpack', 'randomized'}, default='auto') –

  If auto :

  The solver is selected by a default policy based on X.shape and n_components: if the input data is larger than 500x500 and the number of components to extract is lower than 80% of the smallest dimension of the data, then the more efficient ‘randomized’ method is enabled. Otherwise the exact full SVD is computed and optionally truncated afterwards.

  If full :

  run exact full SVD calling the standard LAPACK solver via scipy.linalg.svd and select the components by postprocessing

  If arpack :

  run SVD truncated to n_components calling ARPACK solver via scipy.sparse.linalg.svds. It requires strictly 0 < n_components < min(X.shape)

  If randomized :

  run randomized SVD by the method of Halko et al.

• regressor ({instance of sklearn.kernel_ridge.KernelRidge, precomputed, None}, default=None) –

  The regressor to use for computing the property predictions \(\hat{\mathbf{Y}}\). A pre-fitted regressor may be provided. If the regressor is not None, its kernel parameters (kernel, gamma, degree, coef0, and kernel_params) must be identical to those passed directly to KernelPCovR.

  If precomputed, we assume that the y passed to the fit function is the regressed form of the targets \({\mathbf{\hat{Y}}}\).

• kernel ({'linear', 'poly', 'rbf', 'sigmoid', 'cosine', 'precomputed'} or callable, default='linear') – Kernel.

• gamma (float, default=None) – Kernel coefficient for rbf, poly and sigmoid kernels. Ignored by other kernels.

• degree (int, default=3) – Degree for poly kernels. Ignored by other kernels.

• coef0 (float, default=1) – Independent term in poly and sigmoid kernels. Ignored by other kernels.

• kernel_params (mapping of str to any, default=None) – Parameters (keyword arguments) and values for kernel passed as callable object. Ignored by other kernels.

• center (bool, default=False) – Whether to center any computed kernels

• fit_inverse_transform (bool, default=False) – Learn the inverse transform for non-precomputed kernels. (i.e. learn to find the pre-image of a point)

• tol (float, default=1e-12) – Tolerance for singular values computed by svd_solver == ‘arpack’ and for matrix inversions. Must be of range [0.0, infinity).

• n_jobs (int, default=None) – The number of parallel jobs to run. None means 1 unless in a joblib.parallel_backend context. -1 means using all processors.

• iterated_power (int or 'auto', default='auto') – Number of iterations for the power method computed by svd_solver == ‘randomized’. Must be of range [0, infinity).

• random_state (int, numpy.random.RandomState instance or None, default=None) – Used when the ‘arpack’ or ‘randomized’ solvers are used. Pass an int for reproducible results across multiple function calls.

pt__

pseudo-inverse of the latent-space projection, which can be used to contruct projectors from latent-space

Type:

numpy.darray of size \(({n_{components}, n_{components}})\)

pkt_

the projector, or weights, from the input kernel \(\mathbf{K}\) to the latent-space projection \(\mathbf{T}\)

Type:

numpy.ndarray of size \(({n_{samples}, n_{components}})\)

pky_

the projector, or weights, from the input kernel \(\mathbf{K}\) to the properties \(\mathbf{Y}\)

Type:

numpy.ndarray of size \(({n_{samples}, n_{properties}})\)

pty_

the projector, or weights, from the latent-space projection \(\mathbf{T}\) to the properties \(\mathbf{Y}\)

Type:

numpy.ndarray of size \(({n_{components}, n_{properties}})\)

ptx_

the projector, or weights, from the latent-space projection \(\mathbf{T}\) to the feature matrix \(\mathbf{X}\)

Type:

numpy.ndarray of size \(({n_{components}, n_{features}})\)

X_fit_

The data used to fit the model. This attribute is used to build kernels from new data.

Type:

numpy.ndarray of shape (n_samples, n_features)

Examples

>>> import numpy as np
>>> from skmatter.decomposition import KernelPCovR
>>> from skmatter.preprocessing import StandardFlexibleScaler as SFS
>>> from sklearn.kernel_ridge import KernelRidge
>>> X = np.array([[-1, 1, -3, 1], [1, -2, 1, 2], [-2, 0, -2, -2], [1, 0, 2, -1]])
>>> X = SFS().fit_transform(X)
>>> Y = np.array([[0, -5], [-1, 1], [1, -5], [-3, 2]])
>>> Y = SFS(column_wise=True).fit_transform(Y)
>>> kpcovr = KernelPCovR(
...     mixing=0.1,
...     n_components=2,
...     regressor=KernelRidge(kernel="rbf", gamma=1),
...     kernel="rbf",
...     gamma=1,
... )
>>> kpcovr.fit(X, Y)
KernelPCovR(gamma=1, kernel='rbf', mixing=0.1, n_components=2,
            regressor=KernelRidge(gamma=1, kernel='rbf'))
>>> kpcovr.transform(X)
array([[-0.61261285, -0.18937908],
       [ 0.45242098,  0.25453465],
       [-0.77871824,  0.04847559],
       [ 0.91186937, -0.21211816]])
>>> kpcovr.predict(X)
array([[ 0.5100212 , -0.99488463],
       [-0.18992219,  0.82064368],
       [ 1.11923584, -1.04798016],
       [-1.5635827 ,  1.11078662]])
>>> round(kpcovr.score(X, Y), 5)
np.float64(-0.52039)

fit(X, Y, W=None)

Fit the model with X and Y.

Parameters:

• X (numpy.ndarray, shape (n_samples, n_features)) –

  Training data, where n_samples is the number of samples and n_features is the number of features.

  It is suggested that \(\mathbf{X}\) be centered by its column- means and scaled. If features are related, the matrix should be scaled to have unit variance, otherwise \(\mathbf{X}\) should be scaled so that each feature has a variance of 1 / n_features.

• Y (numpy.ndarray, shape (n_samples, n_properties)) –

  Training data, where n_samples is the number of samples and n_properties is the number of properties

  It is suggested that \(\mathbf{X}\) be centered by its column- means and scaled. If features are related, the matrix should be scaled to have unit variance, otherwise \(\mathbf{Y}\) should be scaled so that each feature has a variance of 1 / n_features.

• W (numpy.ndarray, shape (n_samples, n_properties)) – Regression weights, optional when regressor = precomputed. If not passed, it is assumed that W = np.linalg.lstsq(K, Y, self.tol)[0]

Returns:

self (object) – Returns the instance itself.

transform(X)

Apply dimensionality reduction to X.

X is projected on the first principal components as determined by the modified Kernel PCovR distances.

Parameters:

X (numpy.ndarray, shape (n_samples, n_features)) – New data, where n_samples is the number of samples and n_features is the number of features.

predict(X=None)

Predicts the property values

inverse_transform(T)

Transform input data back to its original space.

\[\mathbf{\hat{X}} = \mathbf{T} \mathbf{P}_{TX} = \mathbf{K} \mathbf{P}_{KT} \mathbf{P}_{TX}\]

Similar to KPCA, the original features are not always recoverable, as the projection is computed from the kernel features, not the original features, and the mapping between the original and kernel features is not one-to-one.

Parameters:

T (numpy.ndarray, shape (n_samples, n_components)) – Projected data, where n_samples is the number of samples and n_components is the number of components.

Returns:

X_original (numpy.ndarray, shape (n_samples, n_features))

score(X, y)

Computes the (negative) loss values for KernelPCovR on the given predictor and response variables.

The loss in \(\mathbf{K}\), as explained in [Helfrecht2020] does not correspond to a traditional Gram loss \(\mathbf{K} - \mathbf{TT}^T\). Indicating the kernel between set A and B as \(\mathbf{K}_{AB}\), the projection of set A as \(\mathbf{T}_A\), and with N and V as the train and validation/test set, one obtains

\[\ell=\frac{\operatorname{Tr}\left[\mathbf{K}_{VV} - 2 \mathbf{K}_{VN} \mathbf{T}_N (\mathbf{T}_N^T \mathbf{T}_N)^{-1} \mathbf{T}_V^T +\mathbf{T}_V(\mathbf{T}_N^T \mathbf{T}_N)^{-1} \mathbf{T}_N^T \mathbf{K}_{NN} \mathbf{T}_N (\mathbf{T}_N^T \mathbf{T}_N)^{-1} \mathbf{T}_V^T\right]}{\operatorname{Tr}(\mathbf{K}_{VV})}\]

The negative loss is returned for easier use in sklearn pipelines, e.g., a grid search, where methods named ‘score’ are meant to be maximized.

Parameters:

• X (numpy.ndarray) – independent (predictor) variable

• Y (numpy.ndarray) – dependent (response) variable

Returns:

L (float) – Negative sum of the KPCA and KRR losses, with the KPCA loss determined by the reconstruction of the kernel

Kernel PCovC

class skmatter.decomposition.KernelPCovC(mixing=0.5, n_components=None, svd_solver='auto', classifier=None, scale_z=False, kernel='linear', gamma=None, degree=3, coef0=1, kernel_params=None, center=False, fit_inverse_transform=False, tol=1e-12, z_mean_tol=1e-12, z_var_tol=1.5, n_jobs=None, iterated_power='auto', random_state=None)

Bases: LinearClassifierMixin, _BaseKPCov

Kernel Principal Covariates Classification (KPCovC).

KPCovC is a modification on the Principal Covariates Classification proposed in [Jorgensen2025]. It determines a latent-space projection \(\mathbf{T}\) which minimizes a combined loss in supervised and unsupervised tasks in the reproducing kernel Hilbert space (RKHS).

This projection is determined by the eigendecomposition of a modified gram matrix \(\mathbf{\tilde{K}}\)

\[\mathbf{\tilde{K}} = \alpha \mathbf{K} + (1 - \alpha) \mathbf{Z}\mathbf{Z}^T\]

where \(\alpha\) is a mixing parameter, \(\mathbf{K}\) is the input kernel of shape \((n_{samples}, n_{samples})\) and \(\mathbf{Z}\) is a matrix of class confidence scores of shape \((n_{samples}, n_{classes})\)

Parameters:

• mixing (float, default=0.5) – mixing parameter, as described in PCovC as \({\alpha}\)

• n_components (int, float or str, default=None) –

  Number of components to keep. if n_components is not set all components are kept:

  n_components == n_samples
  

• svd_solver ({'auto', 'full', 'arpack', 'randomized'}, default='auto') –

  If auto :

  The solver is selected by a default policy based on X.shape and n_components: if the input data is larger than 500x500 and the number of components to extract is lower than 80% of the smallest dimension of the data, then the more efficient ‘randomized’ method is enabled. Otherwise the exact full SVD is computed and optionally truncated afterwards.

  If full :

  run exact full SVD calling the standard LAPACK solver via scipy.linalg.svd and select the components by postprocessing

  If arpack :

  run SVD truncated to n_components calling ARPACK solver via scipy.sparse.linalg.svds. It requires strictly 0 < n_components < min(X.shape)

  If randomized :

  run randomized SVD by the method of Halko et al.

• classifier (estimator object or precomputed, default=None) –

  classifier for computing \({\mathbf{Z}}\). The classifier should be one of the following:

  • sklearn.linear_model.LogisticRegression()

  • sklearn.linear_model.LogisticRegressionCV()

  • sklearn.svm.LinearSVC()

  • sklearn.discriminant_analysis.LinearDiscriminantAnalysis()

  • sklearn.linear_model.RidgeClassifier()

  • sklearn.linear_model.RidgeClassifierCV()

  • sklearn.linear_model.Perceptron()

  If a pre-fitted classifier is provided, it is used to compute \({\mathbf{Z}}\). If None, sklearn.linear_model.LogisticRegression() is used as the classifier.

• scale_z (bool, default=False) – Whether to scale Z prior to eigendecomposition.

• kernel ({"linear", "poly", "rbf", "sigmoid", "precomputed"} or callable, default="linear") – Kernel.

• gamma ({'scale', 'auto'} or float, default=None) – Kernel coefficient for rbf, poly and sigmoid kernels. Ignored by other kernels.

• degree (int, default=3) – Degree for poly kernels. Ignored by other kernels.

• coef0 (float, default=1) – Independent term in poly and sigmoid kernels. Ignored by other kernels.

• kernel_params (mapping of str to any, default=None) – Parameters (keyword arguments) and values for kernel passed as callable object. Ignored by other kernels.

• center (bool, default=False) – Whether to center any computed kernels

• fit_inverse_transform (bool, default=False) – Learn the inverse transform for non-precomputed kernels. (i.e. learn to find the pre-image of a point)

• tol (float, default=1e-12) – Tolerance for singular values computed by svd_solver == ‘arpack’ and for matrix inversions. Must be of range [0.0, infinity).

• z_mean_tol (float, default=1e-12) – Tolerance for the column means of Z. Must be of range [0.0, infinity).

• z_var_tol (float, default=1.5) – Tolerance for the column variances of Z. Must be of range [0.0, infinity).

• n_jobs (int, default=None) – The number of parallel jobs to run. None means 1 unless in a joblib.parallel_backend context. -1 means using all processors.

• iterated_power (int or 'auto', default='auto') – Number of iterations for the power method computed by svd_solver == ‘randomized’. Must be of range [0, infinity).

• random_state (int, numpy.random.RandomState instance or None, default=None) – Used when the ‘arpack’ or ‘randomized’ solvers are used. Pass an int for reproducible results across multiple function calls.

classifier

The linear classifier passed for fitting. If pre-fitted, it is assummed to be fit on a precomputed kernel \(\mathbf{K}\) and \(\mathbf{Y}\).

Type:

estimator object

z_classifier_

The linear classifier fit between the computed kernel \(\mathbf{K}\) and \(\mathbf{Y}\).

Type:

estimator object

classifier_

The linear classifier fit between \(\mathbf{T}\) and \(\mathbf{Y}\).

Type:

estimator object

pt__

pseudo-inverse of the latent-space projection, which can be used to contruct projectors from latent-space

Type:

numpy.darray of size \(({n_{components}, n_{components}})\)

pkt_

the projector, or weights, from the input kernel \(\mathbf{K}\) to the latent-space projection \(\mathbf{T}\)

Type:

numpy.ndarray of size \(({n_{samples}, n_{components}})\)

pkz_

the projector, or weights, from the input kernel \(\mathbf{K}\) to the class confidence scores \(\mathbf{Z}\)

Type:

numpy.ndarray of size \(({n_{samples}, })\) or \(({n_{samples}, n_{classes}})\)

ptz_

the projector, or weights, from the latent-space projection \(\mathbf{T}\) to the class confidence scores \(\mathbf{Z}\)

Type:

numpy.ndarray of size \(({n_{components}, })\) or \(({n_{components}, n_{classes}})\)

ptx_

the projector, or weights, from the latent-space projection \(\mathbf{T}\) to the feature matrix \(\mathbf{X}\)

Type:

numpy.ndarray of size \(({n_{components}, n_{features}})\)

X_fit_

The data used to fit the model. This attribute is used to build kernels from new data.

Type:

numpy.ndarray of shape (n_samples, n_features)

scale_z

Whether Z is being scaled prior to eigendecomposition.

Type:

bool

Examples

>>> import numpy as np
>>> from skmatter.decomposition import KernelPCovC
>>> from sklearn.preprocessing import StandardScaler
>>> X = np.array([[-2, 3, -1, 0], [2, 0, -3, 1], [3, 0, -1, 3], [2, -2, 1, 0]])
>>> X = StandardScaler().fit_transform(X)
>>> Y = np.array([2, 0, 1, 2])
>>> kpcovc = KernelPCovC(
...     mixing=0.1,
...     n_components=2,
...     kernel="rbf",
...     gamma=1,
... )
>>> kpcovc.fit(X, Y)
KernelPCovC(gamma=1, kernel='rbf', mixing=0.1, n_components=2)
>>> kpcovc.transform(X)
array([[-4.45970689e-01,  8.95327566e-06],
       [ 4.52745933e-01,  5.54810948e-01],
       [ 4.52881359e-01, -5.54708315e-01],
       [-4.45921092e-01, -7.32157649e-05]])
>>> kpcovc.predict(X)
array([2, 0, 1, 2])
>>> kpcovc.score(X, Y)
1.0

fit(X, Y, W=None)

Fit the model with X and Y.

A computed kernel K is derived from X, and W is taken from the coefficients of a linear classifier fit between K and Y to compute Z:

\[\mathbf{Z} = \mathbf{K} \mathbf{W}\]

We then call either _fit_feature_space or _fit_sample_space, using Z as our approximation of Y. Finally, we refit a classifier on T and Y to obtain \(\mathbf{P}_{TZ}\).

Parameters:

• X (numpy.ndarray, shape (n_samples, n_features)) –

  Training data, where n_samples is the number of samples and n_features is the number of features.

  It is suggested that \(\mathbf{X}\) be centered by its column- means and scaled. If features are related, the matrix should be scaled to have unit variance, otherwise \(\mathbf{X}\) should be scaled so that each feature has a variance of 1 / n_features.

• Y (numpy.ndarray, shape (n_samples,)) – Training data, where n_samples is the number of samples.

• W (numpy.ndarray, shape (n_features, n_classes)) – Classification weights, optional when classifier = precomputed. If not passed, it is assumed that the weights will be taken from a linear classifier fit between K and Y.

Returns:

self (object) – Returns the instance itself.

transform(X)

Apply dimensionality reduction to X.

X is projected on the first principal components as determined by the modified Kernel PCovR distances.

Parameters:

X (numpy.ndarray, shape (n_samples, n_features)) – New data, where n_samples is the number of samples and n_features is the number of features.

predict(X=None, T=None)

Predicts the property labels using classification on T.

inverse_transform(T)

Transform input data back to its original space.

\[\mathbf{\hat{X}} = \mathbf{T} \mathbf{P}_{TX} = \mathbf{K} \mathbf{P}_{KT} \mathbf{P}_{TX}\]

Similar to KPCA, the original features are not always recoverable, as the projection is computed from the kernel features, not the original features, and the mapping between the original and kernel features is not one-to-one.

Parameters:

T (numpy.ndarray, shape (n_samples, n_components)) – Projected data, where n_samples is the number of samples and n_components is the number of components.

Returns:

X_original (numpy.ndarray, shape (n_samples, n_features))

decision_function(X=None, T=None)

Predicts confidence scores from X or T.

\[\mathbf{Z} = \mathbf{T} \mathbf{P}_{TZ} = \mathbf{K} \mathbf{P}_{KT} \mathbf{P}_{TZ} = \mathbf{K} \mathbf{P}_{KZ}\]

Parameters:

• X (ndarray, shape(n_samples, n_features)) – Original data for which we want to get confidence scores, where n_samples is the number of samples and n_features is the number of features.

• T (ndarray, shape (n_samples, n_components)) – Projected data for which we want to get confidence scores, where n_samples is the number of samples and n_components is the number of components.

Returns:

Z (numpy.ndarray, shape (n_samples,) or (n_samples, n_classes)) – Confidence scores. For binary classification, has shape (n_samples,), for multiclass classification, has shape (n_samples, n_classes)

score(X, y, sample_weight=None)

Return accuracy on provided data and labels.

In multi-label classification, this is the subset accuracy which is a harsh metric since you require for each sample that each label set be correctly predicted.

Parameters:

• X (array-like of shape (n_samples, n_features)) – Test samples.

• y (array-like of shape (n_samples,) or (n_samples, n_outputs)) – True labels for X.

• sample_weight (array-like of shape (n_samples,), default=None) – Sample weights.

Returns:

score (float) – Mean accuracy of self.predict(X) w.r.t. y.

SketchMap

sketch-map [Ceriotti2011] is a nonlinear dimensionality-reduction method for atomistic simulations, where configurations cluster into basins. Short distances mostly measure thermal fluctuations inside a basin, long ones say little more than that two configurations are far apart, so the structure worth preserving lives in between. sketch-map therefore passes the high- and the low-dimensional distances through a saturating sigmoid before matching them, and the embedding is driven by the intermediate range.

The example introduces the method, walks through the landmark workflow used for large datasets, and shows how the implementation reproduces a published sketch-map of the reference C++ code.

class skmatter.decomposition.SketchMap(n_components=2, sigma=None, a_high=None, b_high=None, a_low=None, b_low=None, mds_opt_steps=100, optimizer='L-BFGS-B', max_iter=1000, global_opt_steps='auto', global_optimizer='gradient', mixing_ratio=0.0, mixing_schedule=None, dissimilarity='euclidean', init=None, verbose=False, progress_bar=False)

Bases: TransformerMixin, BaseEstimator

sketch-map embeds high-dimensional data by matching pairwise distances, but only over a chosen range of scales. MDS tries to reproduce every distance and t-SNE keeps only local neighborhoods. sketch-map instead passes both the high- and low-dimensional distances through a sigmoid before comparing them. The sigmoid saturates, so distances far from the switching scale \(\sigma\) contribute little and the fit concentrates on the intermediate range.

Distances much shorter than \(\sigma\) collapse toward 0. In a molecular ensemble these are mostly thermal noise and should not drive the layout. Distances much longer than \(\sigma\) saturate toward 1. In high dimensions long distances are unreliable, so their exact values are dropped while their ordering is kept. Distances near \(\sigma\) pass through almost unchanged and carry the structure the embedding is built around.

The sigmoid is

\[s(r) = 1 - \left(1 + A \left(\frac{r}{\sigma}\right)^a\right)^{-b/a}, \qquad A = 2^{a/b} - 1.\]

Here \(\sigma\) is the switching distance, where \(s(\sigma) = 0.5\). The exponent \(a\) sets how sharply \(s \to 0\) at short range and \(b\) how sharply \(s \to 1\) at long range. The prefactor \(A\) is fixed by \(a\) and \(b\) so the curve always crosses one half at \(\sigma\).

High- and low-dimensional distances use their own exponents (a_high, b_high and a_low, b_low), letting the low dimension take a gentler curve to make room for the structure that cannot fit otherwise.

All five parameters are estimated automatically from the data when left at None. The values actually used are stored in params_ after fitting. If one needs to tune them by hand:

• sigma is the most important parameter, distances below it are treated as “close”, distances above as “far”. Place it just below the peak of the pairwise distance histogram.

• a_high, b_high shape the high-dimensional sigmoid: larger a_high compresses short-range (noise) distances more aggressively, smaller b_high saturates long distances more softly.

• a_low, b_low shape the low-dimensional sigmoid and are usually smaller than their high-dimensional counterparts, to compensate for the volume difference between the spaces.

The fit proceeds through four stages:

1. classical MDS for the initial coordinates (skipped when init is given),

2. refinement of those coordinates against the raw distances (mds_opt_steps),

3. main optimization of the sigmoid-transformed stress (max_iter),

4. global optimization to escape local minima of the non-convex sigmoid stress (global_opt_steps rounds of global_optimizer, with optional mixing annealing via mixing_schedule).

Stages 1-3 always run the same way. The accuracy/cost trade-off lives in stage 4:

• global_optimizer="gradient" (default): repeated L-BFGS relaxations. Fast, any dimensionality.

• global_optimizer="grid": the reference implementation’s pointwise relocation sweeps. Lower stress, 2D only, considerably slower.

• global_optimizer="grid" with mixing_schedule="auto": also anneals the mixing ratio like the reference pipeline. The slowest mode and the one that reproduces published sketch-maps from scratch.

The optimization is deterministic.

Every stage works on the full pairwise distance matrix, so cost and memory grow as \(O(n^2)\). On large datasets the estimator is therefore fit on a set of landmarks, chosen for example by farthest-point sampling and weighted by the population of their Voronoi cells through sample_weight.

Parameters:

• n_components (int, default=2) – Number of dimensions in the target embedding space.

• sigma (float or None, default=None) – Switching distance where \(s(\sigma) = 0.5\), applied to both the high- and low-dimensional distances. If None, it is estimated automatically as 90% of the peak of the pairwise distance distribution. When sample_weight is given the histogram counts each pair with weight \(w_i w_j\), the same statistic the stress sums over, so with Voronoi-weighted landmarks the estimate reflects the full dataset rather than the landmarks uniform spread.

• a_high (float or None, default=None) – Short- and long-range steepness exponents of the high-dimensional sigmoid. If None, both are estimated from the distance distribution.

• b_high (float or None, default=None) – Short- and long-range steepness exponents of the high-dimensional sigmoid. If None, both are estimated from the distance distribution.

• a_low (float or None, default=None) – Short- and long-range steepness exponents of the low-dimensional sigmoid. If None, both are estimated from the input dimensionality.

• b_low (float or None, default=None) – Short- and long-range steepness exponents of the low-dimensional sigmoid. If None, both are estimated from the input dimensionality.

• mds_opt_steps (int, default=100) – Number of stage-2 optimization steps against the raw distances (no sigmoid), refining the classical MDS initialization before the sigmoid stress. Set to 0 to skip the stage. Equivalent to -preopt refinement of the reference C++ implementation.

• optimizer (str, default="L-BFGS-B") – Algorithm used for every full-embedding optimization (stages 2-4). Options: "L-BFGS-B" or "CG".

• max_iter (int, default=1000) – Maximum iterations for the stage-3 main optimization of the (sigmoid-transformed) stress at the fixed mixing_ratio.

• global_opt_steps (int, "auto" or None, default="auto") – Number of global optimization rounds run after the main optimization. The sigmoid stress is non-convex, so the embedding is refined further to escape local minima. What a round does depends on global_optimizer. "auto" (default) uses 5 rounds or up to 10 rounds with early stopping for the annealed global_optimizer="grid". Set to 0 or None to disable.

• global_optimizer ({"gradient", "grid"}, default="gradient") –

  How the global optimization rounds escape local minima.

  "gradient" re-optimizes the whole embedding with L-BFGS, optionally annealing the mixing ratio (see mixing_schedule). It is cheap, works for any n_components, but it only ever moves downhill from the basin it starts in.

  "grid" scans one point at a time over a grid covering the embedding and relocates it when that strictly lowers its stress, briefly optimizing the whole embedding after every accepted move. Moving a single point across the map is a large enough step to cross into another basin, which the gradient rounds cannot do, so it reaches lower stress. It follows the reference C++ implementation. Combined with mixing_schedule="auto" it also anneals the mixing ratio the way the reference C++ pipeline does, which is the closest reproduction of the published sketch-maps this estimator offers. It costs \(O(n^2)\) stress evaluations per sweep, so it is much slower, and it is restricted to n_components=2.

• mixing_ratio (float, default=0.0) –

  Balance between raw distance stress and transformed distance stress:

  • 0.0: Pure sigmoid-transformed stress

  • 1.0: Pure raw distance stress

  • values in between: linear combination

  Used as the constant mixing ratio for the main optimization stage and as the final target of the annealing schedule.

• mixing_schedule (sequence of float, "auto" or None, default=None) –

  Mixing ratios to anneal through during global optimization. None (default) does not anneal and optimizes at the fixed mixing_ratio, matching a single run of the reference C++ implementation, which applies no mixing unless asked.

  With global_optimizer="grid", "auto" replicates the annealing loop of the reference C++ pipeline: the mixing ratio starts at 1 and is multiplied each round by a factor derived from the ratio of the current stresses, so the embedding grows gradually from a small MDS-like map into the sketch-map solution. Because the raw-distance term is typically orders of magnitude larger than the sigmoid term, this error-driven geometric decay is the only way the intermediate mixing values pass through the regime where both terms actually compete.

  With global_optimizer="gradient", each level of the schedule is relaxed with L-BFGS, warm-started from the previous one. "auto" anneals geometrically from 1.0 down to exactly 0 over global_opt_steps levels.

• dissimilarity ({"euclidean", "precomputed"}, default="euclidean") – How X is interpreted in fit(). "euclidean" computes pairwise Euclidean distances from the feature vectors. "precomputed" treats X as a square pairwise distance matrix, so any distance computed elsewhere – periodic, dot-product, or a custom kernel – can be used directly.

• init (array-like of shape (n_samples, n_components) or None, default=None) – Initial embedding coordinates. If None, classical MDS is used.

• verbose (bool, default=False) – If True, print progress information during fitting.

• progress_bar (bool, default=False) – If True, display a tqdm progress bar over the annealing levels of fit(). Requires the optional dependency tqdm.

embedding_

The fitted low-dimensional embedding coordinates.

Type:

ndarray of shape (n_samples, n_components)

stress_

Final stress value (lower is better).

Type:

float

params_

The sigmoid parameters actually used (combination of user-provided and auto-estimated values).

Type:

dict

suggested_params_

The auto-estimated sigmoid parameters, whether or not they were used.

Type:

dict

distance_analysis_

Distance distribution analysis: peak distance, Gaussian range estimates and histogram data.

Type:

dict

n_iter_

Optimizer iterations summed over the full-embedding passes (MDS refinement, main optimization and the per-level relaxations of global optimization).

Type:

int

Examples

Basic usage with automatic parameter estimation:

>>> from skmatter.decomposition import SketchMap
>>> import numpy as np
>>> X = np.random.randn(100, 50)
>>> sm = SketchMap(n_components=2)
>>> embedding = sm.fit_transform(X)
>>> print(embedding.shape)
(100, 2)

Using specific sigmoid parameters:

>>> sm = SketchMap(
...     n_components=2, sigma=7.0, a_high=4.0, b_high=2.0, a_low=2.0, b_low=2.0
... )
>>> embedding = sm.fit_transform(X)
>>> print(embedding.shape)
(100, 2)

Passing a precomputed distance matrix:

>>> from scipy.spatial.distance import cdist
>>> distances = cdist(X, X)
>>> sm = SketchMap(n_components=2, dissimilarity="precomputed")
>>> embedding = sm.fit_transform(distances)
>>> print(embedding.shape)
(100, 2)

Escaping local minima with the grid global optimizer, at a higher cost:

>>> sm = SketchMap(n_components=2, global_optimizer="grid")
>>> embedding = sm.fit_transform(X)
>>> print(embedding.shape)
(100, 2)

Reducing very high-dimensional data with PCA before sketch-map:

>>> from sklearn.decomposition import PCA
>>> from sklearn.pipeline import make_pipeline
>>> pipe = make_pipeline(PCA(n_components=8), SketchMap(n_components=2))
>>> embedding = pipe.fit_transform(X)
>>> print(embedding.shape)
(100, 2)

fit(X, y=None, sample_weight=None)

Fit the sketch-map embedding to the training data

The pairwise distances of X are computed and passed through the sigmoid, the sigmoid parameters are estimated where they were not given and the embedding coordinates are relaxed until the stress is minimized.

Parameters:

• X (array-like of shape (n_samples, n_features), or (n_samples, n_samples)) – Training data. Feature vectors by default, a square pairwise distance matrix when dissimilarity="precomputed".

• y (Ignored) – Not used, present for scikit-learn API compatibility.

• sample_weight (array-like of shape (n_samples,), optional) – Per-sample weights. Samples with higher weights have more influence on the embedding (pair weights are products of sample weights). Default is uniform weights.

Returns:

self (SketchMap) – Returns the fitted instance.

fit_transform(X, y=None, sample_weight=None)

Fit the model and return the embedding.

Parameters:

• X (array-like of shape (n_samples, n_features)) – Training data.

• y (Ignored) – Not used, present for scikit-learn API compatibility.

• sample_weight (array-like of shape (n_samples,), optional) – Per-sample weights.

Returns:

embedding (ndarray of shape (n_samples, n_components)) – Low-dimensional embedding coordinates.