Reference

Setup

ViewDesc

tuple[Hashable, Hashable, Unpack[tuple[Hashable, ...]]]

Main method

solrcmf.SolrCMF

Sparse Orthogonal Low-rank Collective Matrix Factorization (SolrCMF).

Jointly factorises a collection of data matrices X[k], one per view pair k = (i, j, ...), as

X[k] ≈ V[i] diag(d[k]) V[j]^T

where V[i] and V[j] are column-orthonormal factor matrices shared across all data matrices that involve view i or j, and d[k] is a vector of per-factor signal strengths for the view pair k.

Structure sparsity is imposed on d[k] via an L1 penalty (controlled by structure_penalty), causing factors that carry no signal for a given view pair to be zeroed out. Optionally, factor-level sparsity in the loadings can be imposed on a sparse auxiliary variable U via a second L1 penalty (factor_penalty). Globally inactive factors (zero in every d[k]) can be pruned during fitting when factor_pruning=True.

The problem is solved by multi-block ADMM. The algorithm alternates between closed-form updates for Z (data fidelity), D (soft- thresholding), V (Procrustes / SVD), and optionally U and V' (factor sparsity), followed by dual ascent steps for the multipliers.

Attributes:

Name	Type	Description
vs_	dict[Entity, NDArray[float64]]	Fitted orthonormal factor matrices, one per view.
ds_	dict[ViewDesc, NDArray[float64]]	Fitted scaling vectors, one per view pair.
us_	dict[Entity, NDArray[float64]]	Fitted sparse loading matrices, one per view. Only present when factor_penalty is set or factor_pattern is provided.
est_max_rank_	int	Effective rank after fitting (number of globally active factors).

Methods:

Name	Description
__init__	Initialize SolrCMF.
factor_pattern	Return the factor pattern of the fitted solution.
structure_pattern	Return the structure pattern of the fitted solution.

Source code in src/solrcmf/solrcmf.py

class SolrCMF(ADMM[SolrCMFBlocks, SolrCMFConstraints, SolrCMFParams]):
    """Sparse Orthogonal Low-rank Collective Matrix Factorization (SolrCMF).

    Jointly factorises a collection of data matrices X[k], one per view
    pair k = (i, j, ...), as

        X[k] ≈ V[i] diag(d[k]) V[j]^T

    where V[i] and V[j] are column-orthonormal factor matrices shared
    across all data matrices that involve view i or j, and d[k] is a
    vector of per-factor signal strengths for the view pair k.

    Structure sparsity is imposed on d[k] via an L1 penalty (controlled
    by `structure_penalty`), causing factors that carry no signal for a
    given view pair to be zeroed out. Optionally, factor-level sparsity
    in the loadings can be imposed on a sparse auxiliary variable U via a
    second L1 penalty (`factor_penalty`). Globally inactive factors
    (zero in every d[k]) can be pruned during fitting when
    `factor_pruning=True`.

    The problem is solved by multi-block ADMM. The algorithm alternates
    between closed-form updates for Z (data fidelity), D (soft-
    thresholding), V (Procrustes / SVD), and optionally U and V'
    (factor sparsity), followed by dual ascent steps for the multipliers.

    Attributes:
        vs_ (dict[Entity, NDArray[float64]]): Fitted orthonormal factor
            matrices, one per view.
        ds_ (dict[ViewDesc, NDArray[float64]]): Fitted scaling vectors, one per
            view pair.
        us_ (dict[Entity, NDArray[float64]]): Fitted sparse loading matrices,
            one per view. Only present when `factor_penalty` is set or
            `factor_pattern` is provided.
        est_max_rank_ (int): Effective rank after fitting (number of globally
            active factors).

    """

    vs_: dict[Entity, NDArray[float64]]
    ds_: dict[ViewDesc, NDArray[float64]]
    us_: dict[Entity, NDArray[float64]]
    est_max_rank_: int

    _parameter_constraints = {
        **ADMM._parameter_constraints,
        "structure_penalty": [Interval(Real, 0, None, closed="left"), None],
        "max_rank": [Interval(Integral, 1, None, closed="left"), None],
        "factor_penalty": [Interval(Real, 0, None, closed="neither"), None],
        "factor_pruning": ["boolean"],
        "init": [StrOptions({"random", "custom"})],
        "init_kwargs": [dict, None],
        "rho": [Interval(Real, 0.0, None, closed="neither"), None],
        "alpha": [Interval(Real, 0.0, None, closed="left"), None],
        "mu": [Interval(Real, 0.0, None, closed="neither"), None],
    }

    def __init__(
        self,
        *,
        structure_penalty: float | None = None,
        max_rank: int | None = None,
        factor_penalty: float | None = None,
        factor_pruning: bool = True,
        init: str = "random",
        init_kwargs: dict[str, Any] | None = None,
        rho: float | None = None,
        alpha: float | None = None,
        mu: float | None = None,
        max_iter: int = 1000,
        abs_tol: float = 1e-6,
        rel_tol: float = 1e-6,
        save_ctx: bool = False,
    ):
        """Initialize SolrCMF.

        Args:
            structure_penalty: L1 penalty weight on the scaling vectors d[k].
                Larger values produce sparser structure patterns. Mutually
                exclusive with providing `structure_pattern` in `fit`.
            max_rank: Maximum number of factors. Mutually exclusive with
                providing `structure_pattern` in `fit`.
            factor_penalty: L1 penalty weight on the sparse factor loadings U.
                If None, no factor sparsity is imposed. Mutually exclusive with
                providing `factor_pattern` in `fit`.
            factor_pruning: If True, globally inactive factors (zero in all
                d[k]) are removed from all blocks during fitting, reducing the
                effective rank over time.
            init: Initialisation strategy. "random" draws V from the Stiefel
                manifold; "custom" uses factor matrices provided via the `vs`
                and `ds` arguments of `fit`.
            init_kwargs: Additional keyword arguments for the initialiser. For
                "random" init, supports "rng" (int or Generator) for
                reproducibility and "repetitions" (int) for multiple restarts.
            rho: ADMM penalty parameter. If None, a lower bound derived from
                the problem structure is used.
            alpha: Ridge regularisation weight on V. Defaults to 1e-3 * rho.
            mu: Weight for the V' slack penalty. Defaults to 10.0 when factor
                sparsity or a fixed factor pattern is used.
            max_iter: Maximum number of ADMM iterations.
            abs_tol: Absolute convergence tolerance on the objective change.
            rel_tol: Relative convergence tolerance on the objective change.
            save_ctx: If True, the full ADMM context is stored as `ctx_` after
                fitting.

        """
        super().__init__(
            max_iter=max_iter,
            abs_tol=abs_tol,
            rel_tol=rel_tol,
            save_ctx=save_ctx,
        )

        self.structure_penalty = structure_penalty
        self.max_rank = max_rank
        self.factor_penalty = factor_penalty
        self.factor_pruning = factor_pruning
        self.init = init
        self.init_kwargs = init_kwargs
        self.rho = rho
        self.alpha = alpha
        self.mu = mu

    @override
    def _setup(
        self,
        X: dict[ViewDesc, NDArray[float64]],
        *,
        indices: dict[ViewDesc, NDArray[intp]] | None = None,
        structure_weights: (
            dict[ViewDesc, NDArray[float64] | float64] | None
        ) = None,
        structure_pattern: dict[ViewDesc, NDArray[bool_]] | None = None,
        factor_weights: dict[Entity, NDArray[float64] | float64] | None = None,
        factor_pattern: dict[Entity, NDArray[bool_]] | None = None,
        vs: dict[Entity, NDArray[float64]] | None = None,
        ds: dict[ViewDesc, NDArray[float64]] | None = None,
        us: dict[Entity, NDArray[float64]] | None = None,
        **kwargs,
    ):
        # A context will be populated throughout setup
        ctx = Context(SolrCMFBlocks(), SolrCMFConstraints(), SolrCMFParams())

        if not isinstance(X, dict):
            raise TypeError("'X' needs to be a dictionary of data matrices")

        for x in X.values():
            x = check_array(x, ensure_all_finite="allow-nan")

        layout = X.keys()
        # The first two tuple indices indicate the views.
        # The rest are arbitrary to make the indices different if the
        # same views appear.
        views = set([k[i] for k in layout for i in range(2)])
        viewdims_set = set(
            [(k[i], x.shape[i]) for k, x in X.items() for i in range(2)]
        )
        if len(views) != len(viewdims_set):
            raise ValueError(
                "Views do not have consistent dimensions across layout."
                f" Received matrices with dimensions {viewdims_set}."
            )
        viewdims = dict(viewdims_set)

        ctx.data = X

        if not (
            bool(self.structure_penalty is None and self.max_rank is None)
            ^ bool(structure_pattern is None)
        ):
            raise ValueError(
                "Either both structure_penalty and max_rank, or"
                " structure_pattern need to be provided. The respective"
                " other(s) need to be None."
            )

        if self.factor_penalty is not None and factor_pattern is not None:
            raise ValueError(
                "One or both of `factor_penalty` and `factor_pattern`"
                " need to be None."
            )

        if self.structure_penalty is not None and self.max_rank is not None:
            ctx.params.structure_penalty = self.structure_penalty
            max_rank = self.max_rank

            if structure_weights is None:
                ctx.params.structure_weights = {
                    k: float64(1.0) for k in layout
                }
            else:
                # TODO: Add argument check
                ctx.params.structure_weights = structure_weights

            ctx.params.fixed_structure_pattern = False
        else:
            if self.factor_pruning or not isinstance(structure_pattern, dict):
                raise ValueError(
                    "Set 'factor_pruning' to False to use 'structure_pattern'"
                )
            if structure_pattern.keys() != X.keys():
                raise ValueError(
                    "'structure_pattern' must contain one pattern for each"
                    f" data matrix. Expected: {X.keys()}, observed:"
                    f" {structure_pattern.keys()}"
                )

            rks: set[int] = set(p.shape[0] for p in structure_pattern.values())
            if len(rks) != 1:
                raise ValueError(
                    "All patterns in 'structure_pattern' should have the same"
                    f" length. Observed lengths: {rks}"
                )
            # Extract the only element
            max_rank = list(rks)[0]

            ctx.params.structure_pattern = structure_pattern
            ctx.params.fixed_structure_pattern = True

        ctx.params.max_rank = max_rank

        for v, p in viewdims.items():
            if p < max_rank:
                raise ValueError(
                    f"View {v} has dimension {p} which is less than the"
                    f" maximum requested rank {max_rank}"
                )

        if self.factor_penalty is not None:
            if self.mu is None:
                ctx.params.mu = 10.0
            else:
                ctx.params.mu = self.mu

            ctx.params.factor_penalty = self.factor_penalty
            ctx.params.factor_sparsity = True

            if factor_weights is None:
                ctx.params.factor_weights = {
                    k: 1.0 / sqrt(p) for k, p in viewdims.items()
                }
            else:
                # TODO: Add argument check
                ctx.params.factor_weights = factor_weights
        else:
            ctx.params.factor_sparsity = False

        if factor_pattern is not None:
            if self.factor_pruning:
                raise ValueError(
                    "Set 'factor_pruning' to False to use 'factor_pattern'"
                )

            if self.mu is None:
                ctx.params.mu = 10.0
            else:
                ctx.params.mu = self.mu

            # Check factor pattern's correctness
            if factor_pattern.keys() != viewdims.keys():
                raise ValueError(
                    "'factor_pattern' needs to contain a pattern for each"
                    f" view. Views = {viewdims.keys()}, Patterns available"
                    f" for views = {factor_pattern.keys()}"
                )
            dims = {k: p.shape[0] for k, p in factor_pattern.items()}
            if dims != viewdims:
                raise ValueError(
                    f"View dimensions in 'factor_pattern' ({dims}) do not"
                    f" agree with view dimensions in data ({viewdims})."
                )
            rks = set(p.shape[1] for p in factor_pattern.values())
            if len(rks) != 1:
                raise ValueError(
                    "The patterns in 'factor_pattern' need to have the same"
                    f" number of columns. Observed sizes = {rks}"
                )
            # Extract the only element
            rk = list(rks)[0]
            if rk != max_rank:
                raise ValueError(
                    "Number of columns in 'factor_pattern' needs to match"
                    " 'max_rank' or number of elements in each"
                    f" 'structure_pattern'. Expected: {max_rank}, observed:"
                    f" {rk}"
                )

            ctx.params.factor_pattern = factor_pattern
            ctx.params.fixed_factor_pattern = True
        else:
            ctx.params.fixed_factor_pattern = False

        if ctx.params.factor_sparsity or ctx.params.fixed_factor_pattern:
            ctx.params.vp_weights = {
                k: 1.0 / sqrt(p) for k, p in viewdims.items()
            }
            max_vp_w = max([w for w in ctx.params.vp_weights.values()])
            min_vp_w = min([w for w in ctx.params.vp_weights.values()])

            rho_lb = _rho_lower_bound(ctx.params.mu, min_vp_w, max_vp_w)
        else:
            rho_lb = _rho_lower_bound()

        if self.rho is None:
            rho = rho_lb + 0.1
        else:
            rho = self.rho

        if rho <= rho_lb:
            raise ValueError(
                f"rho needs to be greater than {rho_lb} in"
                f" {self.__class__.__name__}; now it is {rho}"
            )

        if ctx.params.factor_sparsity and isinstance(
            self.factor_penalty, float
        ):
            u_edge_cases = {
                k: self.factor_penalty * w * sqrt(viewdims[k]) / rho
                for k, w in ctx.params.factor_weights.items()
            }
            if any([u >= 1 for u in u_edge_cases.values()]):
                warn(
                    "For numerical stability, factor_penalty * weight[k] *"
                    " sqrt(viewdims[k]) / rho < 1 should hold for all views k."
                    f" Here: {u_edge_cases}"
                )

        ctx.params.rho = rho
        if self.alpha is None:
            ctx.params.alpha = 1e-3 * ctx.params.rho
        else:
            ctx.params.alpha = self.alpha

        ctx.params.factor_pruning = self.factor_pruning

        ctx.params.vidx_ridx = {
            v: [(k, k[0]) for k in layout if k[1] == v] for v in views
        }
        ctx.params.vidx_cidx = {
            v: [(k, k[1]) for k in layout if k[0] == v] for v in views
        }

        # Set up ADMM blocks and constraints
        for v in views:
            ctx.add_block("v", v, VBlock, (viewdims[v], max_rank))
        for k in layout:
            ctx.add_block("d", k, DBlock, (max_rank,))

            if self.factor_pruning:
                ctx.blocks.d[k].active_factors = ones((max_rank,), dtype=bool_)

        if self.factor_penalty is not None or factor_pattern is not None:
            for v in views:
                ctx.add_block("u", v, UBlock, (viewdims[v], max_rank))

            # Important to keep vp blocks just before z blocks
            for v in views:
                ctx.add_block("vp", v, VpBlock, (viewdims[v], max_rank))
                ctx.add_constraint(
                    "factor",
                    v,
                    FactorConstraint,
                    (viewdims[v], max_rank),
                )

        # Important to keep z blocks last
        for k in layout:
            ctx.add_block("z", k, ZBlock, (viewdims[k[0]], viewdims[k[1]]))
            ctx.add_constraint(
                "mean_structure",
                k,
                MeanStructureConstraint,
                (viewdims[k[0]], viewdims[k[1]]),
            )

        if self.init_kwargs is None:
            init_kwargs = {}
        else:
            init_kwargs = self.init_kwargs

        # Initialize blocks and constraints
        if self.init == "random":
            init_fn = RandomInitializer(**init_kwargs)
        else:  # self.init == "custom":
            if vs is None or ds is None:
                raise ValueError(
                    f"If 'init' is \"custom\" in {self.__class__.__name__}"
                    " then 'vs' and 'ds' need to be provided to method 'fit'"
                    " as keyword arguments."
                )
            if ctx.params.factor_sparsity or ctx.params.fixed_factor_pattern:
                if us is None:
                    raise ValueError(
                        f"If 'init' is \"custom\" in {self.__class__.__name__}"
                        " then 'us' needs to be provided to method 'fit' as a"
                        " keyword argument."
                    )

            init_fn = FromFormerInitializer(
                vs=vs,
                ds=ds,
                us=us,
                **init_kwargs,
            )

        # Call initializer
        init_fn(ctx)

        # Remove nan entries from `flat_indices` if indices provided.
        # Otherwise have non-nan indices as `flat_indices`
        for k, x in ctx.data.items():
            if indices is None:
                ctx.params.flat_indices[k] = flatnonzero(logical_not(isnan(x)))
            else:
                indices_mask = zeros(x.size, dtype=bool_)
                indices_mask[indices[k]] = True
                not_nan_mask = logical_not(isnan(x)).ravel()
                ctx.params.flat_indices[k] = flatnonzero(
                    logical_and(indices_mask, not_nan_mask)
                )

        return ctx

    @override
    def transform(
        self,
        X: dict[ViewDesc, NDArray[float64]],
        y=None,
        **kwargs,
    ):
        check_is_fitted(self)

        return {
            k: (self.vs_[k[0]] * d) @ self.vs_[k[1]].T
            for k, d in self.ds_.items()
        }

    @override
    def score(
        self,
        X: dict[ViewDesc, NDArray[float64]],
        *,
        indices: dict[ViewDesc, NDArray[intp]] | None = None,
        **kwargs,
    ):
        check_is_fitted(self)

        return neg_mean_squared_error(X, self.transform(X), indices=indices)

    def structure_pattern(self):
        """Return the structure pattern of the fitted solution."""
        check_is_fitted(self)

        return {k: d != 0.0 for k, d in self.ds_.items()}

    def factor_pattern(self) -> dict[Entity, bool] | None:
        """Return the factor pattern of the fitted solution.

        Returns None if factor sparsity was not computed.
        """
        check_is_fitted(self)

        if hasattr(self, "us_"):
            return {k: u != 0.0 for k, u in self.us_.items()}
        else:
            return None

    @override
    def _extra_attrs(
        self, ctx: Context[SolrCMFBlocks, SolrCMFConstraints, SolrCMFParams]
    ):
        out = {}
        out["vs_"] = {k: b.value for k, b in ctx.blocks.v.items()}
        out["ds_"] = {k: b.value for k, b in ctx.blocks.d.items()}
        if ctx.params.factor_sparsity or ctx.params.fixed_factor_pattern:
            out["us_"] = {k: b.value for k, b in ctx.blocks.u.items()}

        out["est_max_rank_"] = sum(
            vstack([d != 0.0 for d in out["ds_"].values()]).sum(0) != 0
        )

        return out

    def _more_tags(self):
        return {
            "X_types": "dict",
        }

__init__(*, structure_penalty=None, max_rank=None, factor_penalty=None, factor_pruning=True, init='random', init_kwargs=None, rho=None, alpha=None, mu=None, max_iter=1000, abs_tol=1e-06, rel_tol=1e-06, save_ctx=False)

Initialize SolrCMF.

Parameters:

Name	Type	Description	Default
structure_penalty	float | None	L1 penalty weight on the scaling vectors d[k]. Larger values produce sparser structure patterns. Mutually exclusive with providing structure_pattern in fit.	None
max_rank	int | None	Maximum number of factors. Mutually exclusive with providing structure_pattern in fit.	None
factor_penalty	float | None	L1 penalty weight on the sparse factor loadings U. If None, no factor sparsity is imposed. Mutually exclusive with providing factor_pattern in fit.	None
factor_pruning	bool	If True, globally inactive factors (zero in all d[k]) are removed from all blocks during fitting, reducing the effective rank over time.	True
init	str	Initialisation strategy. "random" draws V from the Stiefel manifold; "custom" uses factor matrices provided via the vs and ds arguments of fit.	'random'
init_kwargs	dict[str, Any] | None	Additional keyword arguments for the initialiser. For "random" init, supports "rng" (int or Generator) for reproducibility and "repetitions" (int) for multiple restarts.	None
rho	float | None	ADMM penalty parameter. If None, a lower bound derived from the problem structure is used.	None
alpha	float | None	Ridge regularisation weight on V. Defaults to 1e-3 * rho.	None
mu	float | None	Weight for the V' slack penalty. Defaults to 10.0 when factor sparsity or a fixed factor pattern is used.	None
max_iter	int	Maximum number of ADMM iterations.	1000
abs_tol	float	Absolute convergence tolerance on the objective change.	1e-06
rel_tol	float	Relative convergence tolerance on the objective change.	1e-06
save_ctx	bool	If True, the full ADMM context is stored as ctx_ after fitting.	False

Source code in src/solrcmf/solrcmf.py

def __init__(
    self,
    *,
    structure_penalty: float | None = None,
    max_rank: int | None = None,
    factor_penalty: float | None = None,
    factor_pruning: bool = True,
    init: str = "random",
    init_kwargs: dict[str, Any] | None = None,
    rho: float | None = None,
    alpha: float | None = None,
    mu: float | None = None,
    max_iter: int = 1000,
    abs_tol: float = 1e-6,
    rel_tol: float = 1e-6,
    save_ctx: bool = False,
):
    """Initialize SolrCMF.

    Args:
        structure_penalty: L1 penalty weight on the scaling vectors d[k].
            Larger values produce sparser structure patterns. Mutually
            exclusive with providing `structure_pattern` in `fit`.
        max_rank: Maximum number of factors. Mutually exclusive with
            providing `structure_pattern` in `fit`.
        factor_penalty: L1 penalty weight on the sparse factor loadings U.
            If None, no factor sparsity is imposed. Mutually exclusive with
            providing `factor_pattern` in `fit`.
        factor_pruning: If True, globally inactive factors (zero in all
            d[k]) are removed from all blocks during fitting, reducing the
            effective rank over time.
        init: Initialisation strategy. "random" draws V from the Stiefel
            manifold; "custom" uses factor matrices provided via the `vs`
            and `ds` arguments of `fit`.
        init_kwargs: Additional keyword arguments for the initialiser. For
            "random" init, supports "rng" (int or Generator) for
            reproducibility and "repetitions" (int) for multiple restarts.
        rho: ADMM penalty parameter. If None, a lower bound derived from
            the problem structure is used.
        alpha: Ridge regularisation weight on V. Defaults to 1e-3 * rho.
        mu: Weight for the V' slack penalty. Defaults to 10.0 when factor
            sparsity or a fixed factor pattern is used.
        max_iter: Maximum number of ADMM iterations.
        abs_tol: Absolute convergence tolerance on the objective change.
        rel_tol: Relative convergence tolerance on the objective change.
        save_ctx: If True, the full ADMM context is stored as `ctx_` after
            fitting.

    """
    super().__init__(
        max_iter=max_iter,
        abs_tol=abs_tol,
        rel_tol=rel_tol,
        save_ctx=save_ctx,
    )

    self.structure_penalty = structure_penalty
    self.max_rank = max_rank
    self.factor_penalty = factor_penalty
    self.factor_pruning = factor_pruning
    self.init = init
    self.init_kwargs = init_kwargs
    self.rho = rho
    self.alpha = alpha
    self.mu = mu

factor_pattern()

Return the factor pattern of the fitted solution.

Returns None if factor sparsity was not computed.

Source code in src/solrcmf/solrcmf.py

def factor_pattern(self) -> dict[Entity, bool] | None:
    """Return the factor pattern of the fitted solution.

    Returns None if factor sparsity was not computed.
    """
    check_is_fitted(self)

    if hasattr(self, "us_"):
        return {k: u != 0.0 for k, u in self.us_.items()}
    else:
        return None

structure_pattern()

Return the structure pattern of the fitted solution.

Source code in src/solrcmf/solrcmf.py

def structure_pattern(self):
    """Return the structure pattern of the fitted solution."""
    check_is_fitted(self)

    return {k: d != 0.0 for k, d in self.ds_.items()}

Hyperparameter selection

solrcmf.ElementwiseFolds

Element-wise k-fold splitter for collections of data matrices.

Observed entries (non-NaN) across all data matrices are independently partitioned into k roughly equal folds. Each call to split yields k (train_indices, test_indices) pairs, where indices are flat into the corresponding data matrix.

Methods:

Name	Description
__init__	Initialize ElementwiseFolds.
get_n_splits	Return the number of splits.

Source code in src/solrcmf/splits.py

class ElementwiseFolds(BaseSplitter):
    """Element-wise k-fold splitter for collections of data matrices.

    Observed entries (non-NaN) across all data matrices are independently
    partitioned into k roughly equal folds. Each call to `split` yields k
    (train_indices, test_indices) pairs, where indices are flat into the
    corresponding data matrix.
    """

    def __init__(
        self,
        n_splits: int,
        *,
        shuffle: bool = True,
        rng: Generator | None = None,
    ):
        """Initialize ElementwiseFolds.

        Args:
            n_splits: Number of folds. Must be at least 2.
            shuffle: Whether to shuffle observed entries before assigning
                them to folds.
            rng: Random number generator used for shuffling. Must be `None`
                when shuffle is False.

        """
        if n_splits <= 1:
            raise ValueError("n_splits needs to be an integer >= 2")
        self.n_splits = n_splits

        if shuffle is False and rng is not None:
            raise ValueError("rng should be None if shuffle is False")
        self.shuffle = shuffle

        if rng is None:
            rng = default_rng()

        self.rng = rng

    @override
    def _iter_test_indices(self, xs: dict[ViewDesc, NDArray[float64]]):
        """Yield a dict of flat test indices for each fold.

        Observed entries are split into n_splits roughly equal folds
        independently per data matrix.

        """
        # Exclude entries that are already nan
        indices = {
            k: arange(x.size)[flatnonzero(logical_not(isnan(x)))]
            for k, x in xs.items()
        }
        if self.shuffle:
            for idx in indices.values():
                self.rng.shuffle(idx)

        fold_sizes = {
            k: full(self.n_splits, idx.size // self.n_splits, dtype=int_)
            for k, idx in indices.items()
        }
        for k, s in fold_sizes.items():
            s[: indices[k].size % self.n_splits] += 1

        current = {k: 0 for k in fold_sizes.keys()}
        for i in range(self.n_splits):
            test_indices = {
                k: idx[current[k] : current[k] + fold_sizes[k][i]]
                for k, idx in indices.items()
            }
            yield test_indices
            current = {k: idx + fold_sizes[k][i] for k, idx in current.items()}

    @override
    def get_n_splits(self, xs: dict[ViewDesc, NDArray[float64]]):
        """Return the number of splits."""
        return self.n_splits

__init__(n_splits, *, shuffle=True, rng=None)

Initialize ElementwiseFolds.

Parameters:

Name	Type	Description	Default
n_splits	int	Number of folds. Must be at least 2.	required
shuffle	bool	Whether to shuffle observed entries before assigning them to folds.	True
rng	Generator | None	Random number generator used for shuffling. Must be None when shuffle is False.	None

Source code in src/solrcmf/splits.py

def __init__(
    self,
    n_splits: int,
    *,
    shuffle: bool = True,
    rng: Generator | None = None,
):
    """Initialize ElementwiseFolds.

    Args:
        n_splits: Number of folds. Must be at least 2.
        shuffle: Whether to shuffle observed entries before assigning
            them to folds.
        rng: Random number generator used for shuffling. Must be `None`
            when shuffle is False.

    """
    if n_splits <= 1:
        raise ValueError("n_splits needs to be an integer >= 2")
    self.n_splits = n_splits

    if shuffle is False and rng is not None:
        raise ValueError("rng should be None if shuffle is False")
    self.shuffle = shuffle

    if rng is None:
        rng = default_rng()

    self.rng = rng

get_n_splits(xs)

Return the number of splits.

Source code in src/solrcmf/splits.py

@override
def get_n_splits(self, xs: dict[ViewDesc, NDArray[float64]]):
    """Return the number of splits."""
    return self.n_splits

solrcmf.SolrCMFCV

Cross-validated hyperparameter selection for SolrCMF.

Fits SolrCMF over a grid of structure_penalty, max_rank, and factor_penalty values, selects the best combination by cross-validation, and refits on the full data.

Two cross-validation strategies are supported. "structure_first_debiased_cv" first estimates the structure pattern on the full data for each parameter combination, then cross-validates a debiased (unpenalized) refit using that fixed pattern. "penalized_cv" cross-validates the full penalized estimation directly on held-out entries.

Attributes:

Name	Type	Description
cv_results_	dict[str, list]	Per-parameter-combination diagnostics. Each list has one entry per parameter combination. Keys: "structure_penalty", "max_rank", "factor_penalty": the parameter values for each combination. "_fold": held-out score on fold i, where is the chosen scoring function. "mean_", "std_": mean and standard deviation of the per-fold scores across folds. "est_max_rank": effective rank of the best run. "structural_zeros": number of zero entries in d[k] summed over all view pairs (proxy for structure sparsity). "factor_zeros": number of zero entries in U summed over all views (proxy for factor sparsity). Additional keys for cv_strategy="structure_first_debiased_cv": "objective_value_penalized": objective of the penalized fit. "mean_elapsed_process_time_penalized", "std_elapsed_process_time_penalized": CPU time statistics for the penalized structure estimation step. "mean_elapsed_process_time_fixed", "std_elapsed_process_time_fixed": CPU time statistics for the debiased cross-validation step. Additional keys for cv_strategy="penalized_cv": "mean_elapsed_process_time", "std_elapsed_process_time": CPU time statistics across all penalized CV fits.
best_index_	int	Index into cv_results_ of the selected parameter combination.
best_estimator_	SolrCMF	Estimator refit on all data with the selected parameters.
best_max_rank_	int	Effective rank of the best estimator after fitting.

Methods:

Name	Description
__init__	Initialize SolrCMFCV.
fit	Cross-validate the parameter grid and refit on the best combination.

Source code in src/solrcmf/crossval.py

class SolrCMFCV(BaseEstimator):
    """Cross-validated hyperparameter selection for SolrCMF.

    Fits SolrCMF over a grid of `structure_penalty`, `max_rank`, and
    `factor_penalty` values, selects the best combination by
    cross-validation, and refits on the full data.

    Two cross-validation strategies are supported.
    "structure_first_debiased_cv" first estimates the structure pattern on
    the full data for each parameter combination, then cross-validates a
    debiased (unpenalized) refit using that fixed pattern. "penalized_cv"
    cross-validates the full penalized estimation directly on held-out entries.

    Attributes:
        cv_results_ (dict[str, list]): Per-parameter-combination diagnostics.
            Each list has one entry per parameter combination. Keys:

            - "structure_penalty", "max_rank", "factor_penalty":
              the parameter values for each combination.
            - "<score>_fold<i>": held-out score on fold i, where
              <score> is the chosen scoring function.
            - "mean_<score>", "std_<score>": mean and standard
              deviation of the per-fold scores across folds.
            - "est_max_rank": effective rank of the best run.
            - "structural_zeros": number of zero entries in d[k] summed
              over all view pairs (proxy for structure sparsity).
            - "factor_zeros": number of zero entries in U summed over all
              views (proxy for factor sparsity).

            Additional keys for cv_strategy="structure_first_debiased_cv":

            - "objective_value_penalized": objective of the penalized fit.
            - "mean_elapsed_process_time_penalized",
              "std_elapsed_process_time_penalized": CPU time statistics
              for the penalized structure estimation step.
            - "mean_elapsed_process_time_fixed",
              "std_elapsed_process_time_fixed": CPU time statistics for
              the debiased cross-validation step.

            Additional keys for cv_strategy="penalized_cv":

            - "mean_elapsed_process_time",
              "std_elapsed_process_time": CPU time statistics across all
              penalized CV fits.
        best_index_ (int): Index into cv_results_ of the selected parameter
            combination.
        best_estimator_ (SolrCMF): Estimator refit on all data with the
            selected parameters.
        best_max_rank_ (int): Effective rank of the best estimator after
            fitting.

    """

    _parameter_constraints = {
        "structure_penalty": [
            Interval(Real, 0, None, closed="left"),
            "array-like",
        ],
        "max_rank": [Interval(Integral, 1, None, closed="left"), "array-like"],
        "factor_penalty": [
            Interval(Real, 0, None, closed="neither"),
            "array-like",
            None,
        ],
        "factor_pruning": ["boolean"],
        "cv": [Interval(Integral, 2, None, closed="left"), BaseSplitter],
        "cv_strategy": [
            StrOptions({"structure_first_debiased_cv", "penalized_cv"})
        ],
        "score": [
            StrOptions(
                {
                    "neg_mean_squared_error",
                    "neg_sum_squared_error",
                    "weighted_neg_mean_squared_error",
                }
            )
        ],
        "refit": [
            StrOptions(
                {
                    "mean_debiased",
                    "mean_penalized",
                    "1se_debiased",
                    "1se_penalized",
                }
            )
        ],
        "init": [StrOptions({"random", "custom"})],
        "init_kwargs": [dict, None],
        "rho": [Interval(Real, 0.0, None, closed="neither"), None],
        "alpha": [Interval(Real, 0.0, None, closed="left"), None],
        "mu": [Interval(Real, 0.0, None, closed="neither"), None],
        "max_iter": [Interval(Integral, 1, None, closed="left")],
        "abs_tol": [Interval(Real, 0.0, None, closed="neither")],
        "rel_tol": [Interval(Real, 0.0, None, closed="neither")],
        "verbose": ["boolean"],
        "n_jobs": [Integral, None],
    }

    def __init__(
        self,
        *,
        structure_penalty: float | ArrayLike = 1.0,
        max_rank: int | ArrayLike = 10,
        factor_penalty: float | ArrayLike | None = None,
        factor_pruning: bool = True,
        cv: int | BaseSplitter = 10,
        cv_strategy: str = "structure_first_debiased_cv",
        score: str = "neg_mean_squared_error",
        refit: str = "1se_debiased",
        init: str = "random",
        init_kwargs: dict | None = None,
        rho: float | None = None,
        alpha: float | None = None,
        mu: float | None = None,
        max_iter: int = 1000,
        abs_tol: float = 1e-6,
        rel_tol: float = 1e-6,
        verbose: bool = False,
        n_jobs: int | None = None,
    ):
        """Initialize SolrCMFCV.

        Args:
            structure_penalty: L1 penalty weight on d[k], or an array-like of
                values to search over.
            max_rank: Maximum number of factors, or an array-like of values to
                search over.
            factor_penalty: L1 penalty weight on sparse factor loadings U, or
                an array-like of values to search over. `None` disables factor
                sparsity.
            factor_pruning: Whether to prune globally inactive factors during
                fitting.
            cv: Number of cross-validation folds, or a BaseSplitter instance.
            cv_strategy: Cross-validation strategy. See class docstring.
            score: Scoring function used to evaluate held-out predictions.
            refit: Strategy for selecting and refitting the final estimator.
                Prefix "mean" selects by mean CV score; "1se" applies the
                one-standard-error rule. Suffix "debiased" refits without
                penalty using the estimated structure pattern; "penalized"
                refits with the full penalty.
            init: Initialisation strategy passed to each SolrCMF fit.
            init_kwargs: Additional keyword arguments for the initialiser.
            rho: ADMM penalty parameter passed to each SolrCMF fit.
            alpha: Ridge regularisation weight on V passed to each SolrCMF fit.
            mu: V'-slack penalty weight passed to each SolrCMF fit.
            max_iter: Maximum number of ADMM iterations per fit.
            abs_tol: Absolute convergence tolerance per fit.
            rel_tol: Relative convergence tolerance per fit.
            verbose: Whether to print progress messages during fitting.
            n_jobs: Number of parallel jobs. Passed to joblib.Parallel.

        """
        self.structure_penalty = structure_penalty
        self.max_rank = max_rank
        self.factor_penalty = factor_penalty
        self.factor_pruning = factor_pruning
        self.cv = cv
        self.cv_strategy = cv_strategy
        self.score = score
        self.refit = refit
        self.init = init
        self.init_kwargs = init_kwargs
        self.rho = rho
        self.alpha = alpha
        self.mu = mu
        self.max_iter = max_iter
        self.abs_tol = abs_tol
        self.rel_tol = rel_tol
        self.verbose = verbose
        self.n_jobs = n_jobs

    def _check_parameter_grid(self):
        if self.factor_penalty is None:
            factor_penalty = array([None])
        else:
            factor_penalty = atleast_1d(self.factor_penalty)

        # Scalars to 1d-arrays
        structure_penalty, max_rank = atleast_1d(
            self.structure_penalty, self.max_rank
        )

        # Check that all are indeed 1d
        if not (
            ndim(structure_penalty)
            == ndim(max_rank)
            == ndim(factor_penalty)
            == 1
        ):
            raise ValueError(
                f"In {self.__class__.__name__} arguments 'structure_penalty',"
                " 'max_rank', and 'factor_penalty' need to be one-dimensional"
                " or equal to a single number (or 'None' for 'factor_penalty')"
            )

        structure_penalty, max_rank, factor_penalty = broadcast_arrays(
            structure_penalty, max_rank, factor_penalty
        )

        return list(zip(structure_penalty, max_rank, factor_penalty))

    def fit(
        self,
        X: dict[ViewDesc, NDArray[float64]],
        y=None,
        *,
        structure_weights: (
            dict[ViewDesc, NDArray[float64] | float64] | None
        ) = None,
        factor_weights: dict[Entity, NDArray[float64] | float64] | None = None,
        vs: list[dict[Entity, NDArray[float64]]] | None = None,
        ds: list[dict[ViewDesc, NDArray[float64]]] | None = None,
        us: list[dict[Entity, NDArray[float64]]] | None = None,
    ) -> "SolrCMFCV":
        """Cross-validate the parameter grid and refit on the best combination.

        Args:
            X: Data matrices, one per view pair.
            y: Ignored.
            structure_weights: Per-element L1 weights on d[k], one array per
                view pair. Passed through to each SolrCMF fit.
            factor_weights: Per-element L1 weights on U, one array per view.
                Passed through to each SolrCMF fit.
            vs: List of initial factor matrices, one dict per repetition.
                Required when init="custom".
            ds: List of initial scaling vectors, one dict per repetition.
                Required when init="custom".
            us: List of initial sparse loading matrices, one dict per
                repetition. Required when init="custom" and factor sparsity
                is used.

        Returns:
            self

        """
        self._validate_params()

        parameter_grid = self._check_parameter_grid()

        n_params = len(parameter_grid)

        if isinstance(self.cv, Integral):
            cv = ElementwiseFolds(int(self.cv))
        elif isinstance(self.cv, BaseSplitter):
            cv = self.cv

        if self.score == "neg_mean_squared_error":
            score_fn = neg_mean_squared_error
        elif self.score == "neg_sum_squared_error":
            score_fn = neg_sum_squared_error
        elif self.score == "weighted_neg_mean_squared_error":
            score_fn = weighted_neg_mean_squared_error

        results = {
            "structure_penalty": [s for s, _, _ in parameter_grid],
            "max_rank": [m for _, m, _ in parameter_grid],
            "factor_penalty": [f for _, _, f in parameter_grid],
        }

        if self.init_kwargs is None:
            init_kwargs = {}
        else:
            init_kwargs = self.init_kwargs

        if self.init == "random":
            n_reps = 1
            if "repetitions" in init_kwargs:
                n_reps = init_kwargs.pop("repetitions")

            def inits() -> InitsGeneratorType:
                for i in range(n_reps):
                    yield i, (None, None, None)

            # If an rng or seed is supplied, extract it
            if "rng" in init_kwargs:
                rng = default_rng(init_kwargs["rng"])
            else:
                rng = default_rng()
        elif self.init == "custom":
            # If one of these is provided all need to be the same length
            # (if only vs and ds are provided then us is treated as
            # a list of None)
            if not (
                vs is not None and ds is not None and len(vs) == len(ds) >= 1
            ):
                raise ValueError(
                    "If initial values are provided to"
                    f" {self.__class__.__name__}.fit(), then 'vs' and 'ds'"
                    " both need to be provided and have to be the same length"
                )
            if us is not None and len(us) != len(vs):
                raise ValueError(
                    "If initial values for 'u' are provided to"
                    f" {self.__class__.__name__}.fit(), then 'us' needs to"
                    " have the same length as 'vs' and 'ds'"
                )

            n_reps = len(vs)

            def inits() -> InitsGeneratorType:
                for i in range(n_reps):
                    yield i, (vs[i], ds[i], us[i] if us is not None else None)

        else:
            raise ValueError(f"Unknown init method {self.init}")

        base_est = SolrCMF(
            factor_pruning=self.factor_pruning,
            init=self.init,
            init_kwargs=init_kwargs,
            rho=self.rho,
            alpha=self.alpha,
            mu=self.mu,
            max_iter=self.max_iter,
            abs_tol=self.abs_tol,
            rel_tol=self.rel_tol,
        )

        if self.cv_strategy == "structure_first_debiased_cv":
            tmpdir = TemporaryDirectory()
            tmppath = Path(tmpdir.name)

            def _estimate_structure(
                idx_params,
                idx_init,
                structure_penalty,
                max_rank,
                factor_penalty,
                vs,
                ds,
                us,
                rng,
            ):
                est: SolrCMF = clone(base_est)
                est.set_params(
                    structure_penalty=structure_penalty,
                    max_rank=max_rank,
                    factor_penalty=factor_penalty,
                )

                if est.init == "random":
                    est.init_kwargs = {
                        **est.init_kwargs,
                        "rng": default_rng(rng),
                    }

                est.fit(
                    X,
                    structure_weights=structure_weights,
                    factor_weights=factor_weights,
                    vs=vs,
                    ds=ds,
                    us=us,
                )

                if not est.converged_:
                    warn(
                        "Penalized estimation with parameters"
                        f" (structure_penalty={structure_penalty},"
                        f" max_rank={max_rank},"
                        f" factor_penalty={factor_penalty}):"
                        f" {est.__class__.__name__} did not converge"
                        f" after {est.n_iter_} iterations."
                    )

                # Save estimator for later
                dump(est, tmppath / f"{idx_params}_{idx_init}.pkl")

                return (
                    est.objective_value_,
                    est.elapsed_process_time_,
                    est.est_max_rank_,
                    # compute relative to supplied max_rank;
                    # est_max_rank_ could be less in case of
                    # factor pruning
                    (max_rank - est.est_max_rank_) * len(X)
                    + sum(
                        [sum(1 - p) for p in est.structure_pattern().values()]
                    ),
                    (
                        sum(
                            [
                                (max_rank - est.est_max_rank_)
                                * (p.shape[0] - 1)
                                + sum(1 - p)
                                for p in est.factor_pattern().values()
                            ]
                        )
                        if factor_penalty is not None
                        else 0
                    ),
                )

            if self.verbose:
                print(
                    f"Perform structure estimation ({n_reps * n_params} tasks)"
                )

            if self.init == "random":
                # We need to split the randomness for random initialization
                child_states = reshape(
                    rng.bit_generator.spawn(n_reps * n_params),
                    (n_params, n_reps),
                )
            else:
                # Dummy otherwise
                child_states = full((n_params, n_reps), None)

            out = Parallel(
                n_jobs=self.n_jobs, verbose=10 if self.verbose else 0
            )(
                delayed(_estimate_structure)(
                    idx_params,
                    idx_init,
                    structure_penalty,
                    max_rank,
                    factor_penalty,
                    vs,
                    ds,
                    us,
                    child_states[idx_params, idx_init],
                )
                for idx_params, (
                    structure_penalty,
                    max_rank,
                    factor_penalty,
                ) in enumerate(parameter_grid)
                for idx_init, (vs, ds, us) in inits()
            )

            (
                objective_values,
                elapsed_process_times,
                est_max_rank,
                structural_zeros,
                factor_zeros,
            ) = zip(*out)

            if self.verbose:
                print("Determine best runs")

            # Rely on the fact that joblib returns results in the same
            # order as the inputs
            objective_values = split(asarray(objective_values), n_params)
            best_runs = [int(argmin(vals)) for vals in objective_values]
            results["objective_value_penalized"] = [
                vals[idx] for idx, vals in zip(best_runs, objective_values)
            ]

            elapsed_process_times = split(
                asarray(elapsed_process_times), n_params
            )
            results["mean_elapsed_process_time_penalized"] = [
                mean(ts) for ts in elapsed_process_times
            ]
            results["std_elapsed_process_time_penalized"] = [
                std(ts) for ts in elapsed_process_times
            ]

            est_max_rank = split(asarray(est_max_rank), n_params)
            results["est_max_rank"] = [
                rks[idx] for idx, rks in zip(best_runs, est_max_rank)
            ]
            structural_zeros = split(asarray(structural_zeros), n_params)
            results["structural_zeros"] = [
                zs[idx] for idx, zs in zip(best_runs, structural_zeros)
            ]
            factor_zeros = split(asarray(factor_zeros), n_params)
            results["factor_zeros"] = [
                zs[idx] for idx, zs in zip(best_runs, factor_zeros)
            ]

            def _debiased_cv_score(
                est_in: SolrCMF,
                train_indices: dict[ViewDesc, NDArray[intp]],
                test_indices: dict[ViewDesc, NDArray[intp]],
            ):
                est: SolrCMF = clone(base_est)
                est.set_params(
                    init="custom",
                    init_kwargs={"reduce_max_rank": True},
                    factor_pruning=False,  # Set to False always
                )
                est.fit(
                    X,
                    indices=train_indices,
                    structure_pattern=est_in.structure_pattern(),
                    factor_pattern=est_in.factor_pattern(),
                    vs=est_in.vs_,
                    ds=est_in.ds_,
                    us=est_in.us_ if hasattr(est_in, "us_") else None,
                )

                if not est.converged_:
                    warn(
                        "Fixed structure estimation of"
                        f" {est.__class__.__name__} did not converge after"
                        f" {est.n_iter_} iterations."
                    )

                return (
                    score_fn(X, est.transform(X), indices=test_indices),
                    est.elapsed_process_time_,
                )

            # Reads fitted penalized estimators from cache and
            # extracts structure/factor patterns
            def solrcmf_estimators():
                for idx_params, idx_init in zip(range(n_params), best_runs):
                    yield load(tmppath / f"{idx_params}_{idx_init}.pkl")

            # We want exactly the same splits for all parameter combinations,
            # so we produce the splits once and then reuse them.
            cv_splits = list(cv.split(X))
            n_folds = cv.get_n_splits(X)

            if self.verbose:
                print(
                    "Perform debiased cross-validation"
                    f" ({n_params * n_folds} tasks)"
                )

            out = Parallel(
                n_jobs=self.n_jobs, verbose=10 if self.verbose else 0
            )(
                delayed(_debiased_cv_score)(
                    est,
                    train_indices,
                    test_indices,
                )
                for est in solrcmf_estimators()
                for train_indices, test_indices in cv_splits
            )
            (
                scores,
                elapsed_process_times,
            ) = zip(*out)

            for i in range(n_folds):
                results[f"{self.score}_fold{i}"] = [
                    scores[j * n_folds + i] for j in range(n_params)
                ]

            elapsed_process_times = split(
                asarray(elapsed_process_times), n_params
            )
            results["mean_elapsed_process_time_fixed"] = [
                mean(ts) for ts in elapsed_process_times
            ]
            results["std_elapsed_process_time_fixed"] = [
                std(ts) for ts in elapsed_process_times
            ]
        elif self.cv_strategy == "penalized_cv":

            def _penalized_cv_score(
                structure_penalty,
                max_rank,
                factor_penalty,
                vs,
                ds,
                us,
                train_indices,
                test_indices,
                rng,
            ):
                est: SolrCMF = clone(base_est)
                est.set_params(
                    structure_penalty=structure_penalty,
                    max_rank=max_rank,
                    factor_penalty=factor_penalty,
                )

                if est.init == "random":
                    est.init_kwargs = {
                        **est.init_kwargs,
                        "rng": default_rng(rng),
                    }

                est.fit(
                    X,
                    indices=train_indices,
                    structure_weights=structure_weights,
                    factor_weights=factor_weights,
                    vs=vs,
                    ds=ds,
                    us=us,
                )

                if not est.converged_:
                    warn(
                        "Penalized estimation with parameters"
                        f" (structure_penalty={structure_penalty},"
                        f" max_rank={max_rank},"
                        f" factor_penalty={factor_penalty}):"
                        f" {est.__class__.__name__} did not converge"
                        f" after {est.n_iter_} iterations."
                    )

                return (
                    score_fn(X, est.transform(X), indices=test_indices),
                    est.elapsed_process_time_,
                    est.est_max_rank_,
                    # compute relative to supplied max_rank;
                    # est_max_rank_ could be less in case of
                    # factor pruning
                    (max_rank - est.est_max_rank_) * len(X)
                    + sum(
                        [sum(1 - p) for p in est.structure_pattern().values()]
                    ),
                    (
                        sum(
                            [
                                (max_rank - est.est_max_rank_)
                                * (p.shape[0] - 1)
                                + sum(1 - p)
                                for p in est.factor_pattern().values()
                            ]
                        )
                        if factor_penalty is not None
                        else 0
                    ),
                )

            # We want exactly the same splits for all parameter combinations,
            # so we produce the splits once and then reuse them.
            cv_splits = list(cv.split(X))
            n_folds = cv.get_n_splits(X)

            if self.verbose:
                print(
                    "Perform penalized cross-validation"
                    f" ({n_params * n_reps * n_folds} tasks)"
                )

            if self.init == "random":
                # We need to split the randomness for random initialization
                child_states = reshape(
                    rng.bit_generator.spawn(n_reps * n_params * n_folds),
                    (n_params, n_reps, n_folds),
                )
            else:
                # Dummy otherwise
                child_states = full((n_params, n_reps, n_folds), None)

            out = Parallel(
                n_jobs=self.n_jobs, verbose=10 if self.verbose else 0
            )(
                delayed(_penalized_cv_score)(
                    structure_penalty,
                    max_rank,
                    factor_penalty,
                    vs,
                    ds,
                    us,
                    train_indices,
                    test_indices,
                    child_states[idx_param, idx_init, idx_fold],
                )
                for idx_param, (
                    structure_penalty,
                    max_rank,
                    factor_penalty,
                ) in enumerate(parameter_grid)
                for idx_init, (vs, ds, us) in inits()
                for idx_fold, (train_indices, test_indices) in enumerate(
                    cv_splits
                )
            )

            (
                scores,
                elapsed_process_times,
                est_max_rank,
                structural_zeros,
                factor_zeros,
            ) = zip(*out)

            for i in range(n_folds):
                results[f"{self.score}_fold{i}"] = [nan] * n_params

            best_runs = [-1] * n_params
            best_score = [inf] * n_params
            for idx_params, scores_params in enumerate(
                split(asarray(scores), n_params)
            ):
                for idx_init, scores_inits in enumerate(
                    split(scores_params, n_reps)
                ):
                    if mean(scores_inits) < best_score[idx_params]:
                        best_score[idx_params] = mean(scores_inits)
                        best_runs[idx_params] = idx_init
                        for i in range(n_folds):
                            results[f"{self.score}_fold{i}"][idx_params] = (
                                scores_inits[i]
                            )

            elapsed_process_times = split(
                asarray(elapsed_process_times), n_params
            )
            results["mean_elapsed_process_time"] = [
                mean(ts) for ts in elapsed_process_times
            ]
            results["std_elapsed_process_time"] = [
                std(ts) for ts in elapsed_process_times
            ]

            results["est_max_rank"] = [
                mean(
                    est_max_rank[
                        (
                            idx_params * n_reps * n_folds
                            + best_runs[idx_params] * n_folds
                        ) : (
                            idx_params * n_reps * n_folds
                            + (best_runs[idx_params] + 1) * n_folds
                        )
                    ]
                )
                for idx_params in range(n_params)
            ]
            results["structural_zeros"] = [
                mean(
                    structural_zeros[
                        (
                            idx_params * n_reps * n_folds
                            + best_runs[idx_params] * n_folds
                        ) : (
                            idx_params * n_reps * n_folds
                            + (best_runs[idx_params] + 1) * n_folds
                        )
                    ]
                )
                for idx_params in range(n_params)
            ]

            results["factor_zeros"] = [
                mean(
                    factor_zeros[
                        (
                            idx_params * n_reps * n_folds
                            + best_runs[idx_params] * n_folds
                        ) : (
                            idx_params * n_reps * n_folds
                            + (best_runs[idx_params] + 1) * n_folds
                        )
                    ]
                )
                for idx_params in range(n_params)
            ]

        # Post-processing on the full dictionary. Same for both cases
        scores = vstack(
            [results[f"{self.score}_fold{i}"] for i in range(n_folds)]
        )
        results.update(
            {
                f"mean_{self.score}": scores.mean(0),
                f"std_{self.score}": scores.std(0),
            }
        )

        self.cv_results_ = results

        if self.verbose:
            print("Re-fit final estimator")

        if self.refit.startswith("mean"):
            self.best_index_ = argmax(results[f"mean_{self.score}"])
        elif self.refit.startswith("1se"):
            # Choose the solution with maximal structure sparsity within
            # 1 standard error of the best solution
            max_index = argmax(results[f"mean_{self.score}"])

            candidates = flatnonzero(
                results[f"mean_{self.score}"]
                >= (
                    results[f"mean_{self.score}"][max_index]
                    - results[f"std_{self.score}"][max_index]
                )
            )

            # Primarily choose the solution with the most
            # structural zeros and then select the solution with the
            # most factor zeros if factor sparsity was requested
            structural_zeros = asarray(
                [results["structural_zeros"][i] for i in candidates]
            )
            most_sz_candidates = candidates[
                flatnonzero(structural_zeros == max(structural_zeros))
            ]

            factor_zeros = [
                results["factor_zeros"][i] for i in most_sz_candidates
            ]

            self.best_index_ = most_sz_candidates[argmax(factor_zeros)]

        structure_penalty, max_rank, factor_penalty = parameter_grid[
            self.best_index_
        ]

        if self.verbose:
            print(
                "Best fit with\n"
                f"  structure_penalty = {structure_penalty}\n"
                f"  max_rank = {max_rank}\n"
                f"  factor_penalty = {factor_penalty}\n\n"
                "  estimated max_rank = "
                f"{results['est_max_rank'][self.best_index_]}"
            )

        # Re-fit best run on all data
        if self.cv_strategy == "structure_first_debiased_cv":
            # Load respective penalized estimator from cache
            est = load(
                tmppath
                / f"{self.best_index_}_{best_runs[self.best_index_]}.pkl"
            )

            self.best_estimator_: SolrCMF = clone(base_est)

            if self.refit.endswith("debiased"):
                self.best_estimator_.set_params(
                    init="custom",
                    init_kwargs={"reduce_max_rank": True},
                    factor_pruning=False,  # Set to False always
                )
                self.best_estimator_.fit(
                    X,
                    structure_pattern=est.structure_pattern(),
                    factor_pattern=est.factor_pattern(),
                    vs=est.vs_,
                    ds=est.ds_,
                    us=est.us_ if hasattr(est, "us_") else None,
                )
            elif self.refit.endswith("penalized"):
                self.best_estimator_.set_params(
                    structure_penalty=structure_penalty,
                    max_rank=max_rank,
                    factor_penalty=factor_penalty,
                    init="custom",
                )
                self.best_estimator_.fit(
                    X,
                    vs=est.vs_,
                    ds=est.ds_,
                    us=est.us_ if hasattr(est, "us_") else None,
                )

            tmpdir.cleanup()
        elif self.cv_strategy == "penalized_cv":
            # A penalized fit needs to be performed irrespectively.
            # Either because it is the final fit or because we need the
            # structure/factor pattern.
            if self.init == "custom":
                assert vs is not None and ds is not None
                vs_init = vs[best_runs[self.best_index_]]
                ds_init = ds[best_runs[self.best_index_]]
                if us is not None:
                    us_init = us[best_runs[self.best_index_]]
                else:
                    us_init = None
            else:
                vs_init = None
                ds_init = None
                us_init = None

            final_est: SolrCMF = clone(base_est)
            final_est.set_params(
                structure_penalty=structure_penalty,
                max_rank=max_rank,
                factor_penalty=factor_penalty,
            )
            final_est.fit(
                X,
                vs=vs_init,
                ds=ds_init,
                us=us_init,
            )

            if self.refit.endswith("debiased"):
                final_est_debiased: SolrCMF = clone(base_est)
                final_est_debiased.set_params(
                    init="custom",
                    init_kwargs={"reduce_max_rank": True},
                    factor_pruning=False,  # Set to False always
                )
                final_est_debiased.fit(
                    X,
                    structure_pattern=final_est.structure_pattern(),
                    factor_pattern=final_est.factor_pattern(),
                    vs=final_est.vs_,
                    ds=final_est.ds_,
                    us=final_est.us_ if hasattr(final_est, "us_") else None,
                )

                self.best_estimator_ = final_est_debiased
            elif self.refit.endswith("penalized"):
                self.best_estimator_ = final_est

        self.best_max_rank_ = self.best_estimator_.est_max_rank_

        return self

__init__(*, structure_penalty=1.0, max_rank=10, factor_penalty=None, factor_pruning=True, cv=10, cv_strategy='structure_first_debiased_cv', score='neg_mean_squared_error', refit='1se_debiased', init='random', init_kwargs=None, rho=None, alpha=None, mu=None, max_iter=1000, abs_tol=1e-06, rel_tol=1e-06, verbose=False, n_jobs=None)

Initialize SolrCMFCV.

Parameters:

Name	Type	Description	Default
structure_penalty	float | ArrayLike	L1 penalty weight on d[k], or an array-like of values to search over.	1.0
max_rank	int | ArrayLike	Maximum number of factors, or an array-like of values to search over.	10
factor_penalty	float | ArrayLike | None	L1 penalty weight on sparse factor loadings U, or an array-like of values to search over. None disables factor sparsity.	None
factor_pruning	bool	Whether to prune globally inactive factors during fitting.	True
cv	int | BaseSplitter	Number of cross-validation folds, or a BaseSplitter instance.	10
cv_strategy	str	Cross-validation strategy. See class docstring.	'structure_first_debiased_cv'
score	str	Scoring function used to evaluate held-out predictions.	'neg_mean_squared_error'
refit	str	Strategy for selecting and refitting the final estimator. Prefix "mean" selects by mean CV score; "1se" applies the one-standard-error rule. Suffix "debiased" refits without penalty using the estimated structure pattern; "penalized" refits with the full penalty.	'1se_debiased'
init	str	Initialisation strategy passed to each SolrCMF fit.	'random'
init_kwargs	dict | None	Additional keyword arguments for the initialiser.	None
rho	float | None	ADMM penalty parameter passed to each SolrCMF fit.	None
alpha	float | None	Ridge regularisation weight on V passed to each SolrCMF fit.	None
mu	float | None	V'-slack penalty weight passed to each SolrCMF fit.	None
max_iter	int	Maximum number of ADMM iterations per fit.	1000
abs_tol	float	Absolute convergence tolerance per fit.	1e-06
rel_tol	float	Relative convergence tolerance per fit.	1e-06
verbose	bool	Whether to print progress messages during fitting.	False
n_jobs	int | None	Number of parallel jobs. Passed to joblib.Parallel.	None

Source code in src/solrcmf/crossval.py

def __init__(
    self,
    *,
    structure_penalty: float | ArrayLike = 1.0,
    max_rank: int | ArrayLike = 10,
    factor_penalty: float | ArrayLike | None = None,
    factor_pruning: bool = True,
    cv: int | BaseSplitter = 10,
    cv_strategy: str = "structure_first_debiased_cv",
    score: str = "neg_mean_squared_error",
    refit: str = "1se_debiased",
    init: str = "random",
    init_kwargs: dict | None = None,
    rho: float | None = None,
    alpha: float | None = None,
    mu: float | None = None,
    max_iter: int = 1000,
    abs_tol: float = 1e-6,
    rel_tol: float = 1e-6,
    verbose: bool = False,
    n_jobs: int | None = None,
):
    """Initialize SolrCMFCV.

    Args:
        structure_penalty: L1 penalty weight on d[k], or an array-like of
            values to search over.
        max_rank: Maximum number of factors, or an array-like of values to
            search over.
        factor_penalty: L1 penalty weight on sparse factor loadings U, or
            an array-like of values to search over. `None` disables factor
            sparsity.
        factor_pruning: Whether to prune globally inactive factors during
            fitting.
        cv: Number of cross-validation folds, or a BaseSplitter instance.
        cv_strategy: Cross-validation strategy. See class docstring.
        score: Scoring function used to evaluate held-out predictions.
        refit: Strategy for selecting and refitting the final estimator.
            Prefix "mean" selects by mean CV score; "1se" applies the
            one-standard-error rule. Suffix "debiased" refits without
            penalty using the estimated structure pattern; "penalized"
            refits with the full penalty.
        init: Initialisation strategy passed to each SolrCMF fit.
        init_kwargs: Additional keyword arguments for the initialiser.
        rho: ADMM penalty parameter passed to each SolrCMF fit.
        alpha: Ridge regularisation weight on V passed to each SolrCMF fit.
        mu: V'-slack penalty weight passed to each SolrCMF fit.
        max_iter: Maximum number of ADMM iterations per fit.
        abs_tol: Absolute convergence tolerance per fit.
        rel_tol: Relative convergence tolerance per fit.
        verbose: Whether to print progress messages during fitting.
        n_jobs: Number of parallel jobs. Passed to joblib.Parallel.

    """
    self.structure_penalty = structure_penalty
    self.max_rank = max_rank
    self.factor_penalty = factor_penalty
    self.factor_pruning = factor_pruning
    self.cv = cv
    self.cv_strategy = cv_strategy
    self.score = score
    self.refit = refit
    self.init = init
    self.init_kwargs = init_kwargs
    self.rho = rho
    self.alpha = alpha
    self.mu = mu
    self.max_iter = max_iter
    self.abs_tol = abs_tol
    self.rel_tol = rel_tol
    self.verbose = verbose
    self.n_jobs = n_jobs

fit(X, y=None, *, structure_weights=None, factor_weights=None, vs=None, ds=None, us=None)

Cross-validate the parameter grid and refit on the best combination.

Parameters:

Name	Type	Description	Default
X	dict[ViewDesc, NDArray[float64]]	Data matrices, one per view pair.	required
y		Ignored.	None
structure_weights	dict[ViewDesc, NDArray[float64] | float64] | None	Per-element L1 weights on d[k], one array per view pair. Passed through to each SolrCMF fit.	None
factor_weights	dict[Entity, NDArray[float64] | float64] | None	Per-element L1 weights on U, one array per view. Passed through to each SolrCMF fit.	None
vs	list[dict[Entity, NDArray[float64]]] | None	List of initial factor matrices, one dict per repetition. Required when init="custom".	None
ds	list[dict[ViewDesc, NDArray[float64]]] | None	List of initial scaling vectors, one dict per repetition. Required when init="custom".	None
us	list[dict[Entity, NDArray[float64]]] | None	List of initial sparse loading matrices, one dict per repetition. Required when init="custom" and factor sparsity is used.	None

Returns:

Type	Description
SolrCMFCV	self

Source code in src/solrcmf/crossval.py

def fit(
    self,
    X: dict[ViewDesc, NDArray[float64]],
    y=None,
    *,
    structure_weights: (
        dict[ViewDesc, NDArray[float64] | float64] | None
    ) = None,
    factor_weights: dict[Entity, NDArray[float64] | float64] | None = None,
    vs: list[dict[Entity, NDArray[float64]]] | None = None,
    ds: list[dict[ViewDesc, NDArray[float64]]] | None = None,
    us: list[dict[Entity, NDArray[float64]]] | None = None,
) -> "SolrCMFCV":
    """Cross-validate the parameter grid and refit on the best combination.

    Args:
        X: Data matrices, one per view pair.
        y: Ignored.
        structure_weights: Per-element L1 weights on d[k], one array per
            view pair. Passed through to each SolrCMF fit.
        factor_weights: Per-element L1 weights on U, one array per view.
            Passed through to each SolrCMF fit.
        vs: List of initial factor matrices, one dict per repetition.
            Required when init="custom".
        ds: List of initial scaling vectors, one dict per repetition.
            Required when init="custom".
        us: List of initial sparse loading matrices, one dict per
            repetition. Required when init="custom" and factor sparsity
            is used.

    Returns:
        self

    """
    self._validate_params()

    parameter_grid = self._check_parameter_grid()

    n_params = len(parameter_grid)

    if isinstance(self.cv, Integral):
        cv = ElementwiseFolds(int(self.cv))
    elif isinstance(self.cv, BaseSplitter):
        cv = self.cv

    if self.score == "neg_mean_squared_error":
        score_fn = neg_mean_squared_error
    elif self.score == "neg_sum_squared_error":
        score_fn = neg_sum_squared_error
    elif self.score == "weighted_neg_mean_squared_error":
        score_fn = weighted_neg_mean_squared_error

    results = {
        "structure_penalty": [s for s, _, _ in parameter_grid],
        "max_rank": [m for _, m, _ in parameter_grid],
        "factor_penalty": [f for _, _, f in parameter_grid],
    }

    if self.init_kwargs is None:
        init_kwargs = {}
    else:
        init_kwargs = self.init_kwargs

    if self.init == "random":
        n_reps = 1
        if "repetitions" in init_kwargs:
            n_reps = init_kwargs.pop("repetitions")

        def inits() -> InitsGeneratorType:
            for i in range(n_reps):
                yield i, (None, None, None)

        # If an rng or seed is supplied, extract it
        if "rng" in init_kwargs:
            rng = default_rng(init_kwargs["rng"])
        else:
            rng = default_rng()
    elif self.init == "custom":
        # If one of these is provided all need to be the same length
        # (if only vs and ds are provided then us is treated as
        # a list of None)
        if not (
            vs is not None and ds is not None and len(vs) == len(ds) >= 1
        ):
            raise ValueError(
                "If initial values are provided to"
                f" {self.__class__.__name__}.fit(), then 'vs' and 'ds'"
                " both need to be provided and have to be the same length"
            )
        if us is not None and len(us) != len(vs):
            raise ValueError(
                "If initial values for 'u' are provided to"
                f" {self.__class__.__name__}.fit(), then 'us' needs to"
                " have the same length as 'vs' and 'ds'"
            )

        n_reps = len(vs)

        def inits() -> InitsGeneratorType:
            for i in range(n_reps):
                yield i, (vs[i], ds[i], us[i] if us is not None else None)

    else:
        raise ValueError(f"Unknown init method {self.init}")

    base_est = SolrCMF(
        factor_pruning=self.factor_pruning,
        init=self.init,
        init_kwargs=init_kwargs,
        rho=self.rho,
        alpha=self.alpha,
        mu=self.mu,
        max_iter=self.max_iter,
        abs_tol=self.abs_tol,
        rel_tol=self.rel_tol,
    )

    if self.cv_strategy == "structure_first_debiased_cv":
        tmpdir = TemporaryDirectory()
        tmppath = Path(tmpdir.name)

        def _estimate_structure(
            idx_params,
            idx_init,
            structure_penalty,
            max_rank,
            factor_penalty,
            vs,
            ds,
            us,
            rng,
        ):
            est: SolrCMF = clone(base_est)
            est.set_params(
                structure_penalty=structure_penalty,
                max_rank=max_rank,
                factor_penalty=factor_penalty,
            )

            if est.init == "random":
                est.init_kwargs = {
                    **est.init_kwargs,
                    "rng": default_rng(rng),
                }

            est.fit(
                X,
                structure_weights=structure_weights,
                factor_weights=factor_weights,
                vs=vs,
                ds=ds,
                us=us,
            )

            if not est.converged_:
                warn(
                    "Penalized estimation with parameters"
                    f" (structure_penalty={structure_penalty},"
                    f" max_rank={max_rank},"
                    f" factor_penalty={factor_penalty}):"
                    f" {est.__class__.__name__} did not converge"
                    f" after {est.n_iter_} iterations."
                )

            # Save estimator for later
            dump(est, tmppath / f"{idx_params}_{idx_init}.pkl")

            return (
                est.objective_value_,
                est.elapsed_process_time_,
                est.est_max_rank_,
                # compute relative to supplied max_rank;
                # est_max_rank_ could be less in case of
                # factor pruning
                (max_rank - est.est_max_rank_) * len(X)
                + sum(
                    [sum(1 - p) for p in est.structure_pattern().values()]
                ),
                (
                    sum(
                        [
                            (max_rank - est.est_max_rank_)
                            * (p.shape[0] - 1)
                            + sum(1 - p)
                            for p in est.factor_pattern().values()
                        ]
                    )
                    if factor_penalty is not None
                    else 0
                ),
            )

        if self.verbose:
            print(
                f"Perform structure estimation ({n_reps * n_params} tasks)"
            )

        if self.init == "random":
            # We need to split the randomness for random initialization
            child_states = reshape(
                rng.bit_generator.spawn(n_reps * n_params),
                (n_params, n_reps),
            )
        else:
            # Dummy otherwise
            child_states = full((n_params, n_reps), None)

        out = Parallel(
            n_jobs=self.n_jobs, verbose=10 if self.verbose else 0
        )(
            delayed(_estimate_structure)(
                idx_params,
                idx_init,
                structure_penalty,
                max_rank,
                factor_penalty,
                vs,
                ds,
                us,
                child_states[idx_params, idx_init],
            )
            for idx_params, (
                structure_penalty,
                max_rank,
                factor_penalty,
            ) in enumerate(parameter_grid)
            for idx_init, (vs, ds, us) in inits()
        )

        (
            objective_values,
            elapsed_process_times,
            est_max_rank,
            structural_zeros,
            factor_zeros,
        ) = zip(*out)

        if self.verbose:
            print("Determine best runs")

        # Rely on the fact that joblib returns results in the same
        # order as the inputs
        objective_values = split(asarray(objective_values), n_params)
        best_runs = [int(argmin(vals)) for vals in objective_values]
        results["objective_value_penalized"] = [
            vals[idx] for idx, vals in zip(best_runs, objective_values)
        ]

        elapsed_process_times = split(
            asarray(elapsed_process_times), n_params
        )
        results["mean_elapsed_process_time_penalized"] = [
            mean(ts) for ts in elapsed_process_times
        ]
        results["std_elapsed_process_time_penalized"] = [
            std(ts) for ts in elapsed_process_times
        ]

        est_max_rank = split(asarray(est_max_rank), n_params)
        results["est_max_rank"] = [
            rks[idx] for idx, rks in zip(best_runs, est_max_rank)
        ]
        structural_zeros = split(asarray(structural_zeros), n_params)
        results["structural_zeros"] = [
            zs[idx] for idx, zs in zip(best_runs, structural_zeros)
        ]
        factor_zeros = split(asarray(factor_zeros), n_params)
        results["factor_zeros"] = [
            zs[idx] for idx, zs in zip(best_runs, factor_zeros)
        ]

        def _debiased_cv_score(
            est_in: SolrCMF,
            train_indices: dict[ViewDesc, NDArray[intp]],
            test_indices: dict[ViewDesc, NDArray[intp]],
        ):
            est: SolrCMF = clone(base_est)
            est.set_params(
                init="custom",
                init_kwargs={"reduce_max_rank": True},
                factor_pruning=False,  # Set to False always
            )
            est.fit(
                X,
                indices=train_indices,
                structure_pattern=est_in.structure_pattern(),
                factor_pattern=est_in.factor_pattern(),
                vs=est_in.vs_,
                ds=est_in.ds_,
                us=est_in.us_ if hasattr(est_in, "us_") else None,
            )

            if not est.converged_:
                warn(
                    "Fixed structure estimation of"
                    f" {est.__class__.__name__} did not converge after"
                    f" {est.n_iter_} iterations."
                )

            return (
                score_fn(X, est.transform(X), indices=test_indices),
                est.elapsed_process_time_,
            )

        # Reads fitted penalized estimators from cache and
        # extracts structure/factor patterns
        def solrcmf_estimators():
            for idx_params, idx_init in zip(range(n_params), best_runs):
                yield load(tmppath / f"{idx_params}_{idx_init}.pkl")

        # We want exactly the same splits for all parameter combinations,
        # so we produce the splits once and then reuse them.
        cv_splits = list(cv.split(X))
        n_folds = cv.get_n_splits(X)

        if self.verbose:
            print(
                "Perform debiased cross-validation"
                f" ({n_params * n_folds} tasks)"
            )

        out = Parallel(
            n_jobs=self.n_jobs, verbose=10 if self.verbose else 0
        )(
            delayed(_debiased_cv_score)(
                est,
                train_indices,
                test_indices,
            )
            for est in solrcmf_estimators()
            for train_indices, test_indices in cv_splits
        )
        (
            scores,
            elapsed_process_times,
        ) = zip(*out)

        for i in range(n_folds):
            results[f"{self.score}_fold{i}"] = [
                scores[j * n_folds + i] for j in range(n_params)
            ]

        elapsed_process_times = split(
            asarray(elapsed_process_times), n_params
        )
        results["mean_elapsed_process_time_fixed"] = [
            mean(ts) for ts in elapsed_process_times
        ]
        results["std_elapsed_process_time_fixed"] = [
            std(ts) for ts in elapsed_process_times
        ]
    elif self.cv_strategy == "penalized_cv":

        def _penalized_cv_score(
            structure_penalty,
            max_rank,
            factor_penalty,
            vs,
            ds,
            us,
            train_indices,
            test_indices,
            rng,
        ):
            est: SolrCMF = clone(base_est)
            est.set_params(
                structure_penalty=structure_penalty,
                max_rank=max_rank,
                factor_penalty=factor_penalty,
            )

            if est.init == "random":
                est.init_kwargs = {
                    **est.init_kwargs,
                    "rng": default_rng(rng),
                }

            est.fit(
                X,
                indices=train_indices,
                structure_weights=structure_weights,
                factor_weights=factor_weights,
                vs=vs,
                ds=ds,
                us=us,
            )

            if not est.converged_:
                warn(
                    "Penalized estimation with parameters"
                    f" (structure_penalty={structure_penalty},"
                    f" max_rank={max_rank},"
                    f" factor_penalty={factor_penalty}):"
                    f" {est.__class__.__name__} did not converge"
                    f" after {est.n_iter_} iterations."
                )

            return (
                score_fn(X, est.transform(X), indices=test_indices),
                est.elapsed_process_time_,
                est.est_max_rank_,
                # compute relative to supplied max_rank;
                # est_max_rank_ could be less in case of
                # factor pruning
                (max_rank - est.est_max_rank_) * len(X)
                + sum(
                    [sum(1 - p) for p in est.structure_pattern().values()]
                ),
                (
                    sum(
                        [
                            (max_rank - est.est_max_rank_)
                            * (p.shape[0] - 1)
                            + sum(1 - p)
                            for p in est.factor_pattern().values()
                        ]
                    )
                    if factor_penalty is not None
                    else 0
                ),
            )

        # We want exactly the same splits for all parameter combinations,
        # so we produce the splits once and then reuse them.
        cv_splits = list(cv.split(X))
        n_folds = cv.get_n_splits(X)

        if self.verbose:
            print(
                "Perform penalized cross-validation"
                f" ({n_params * n_reps * n_folds} tasks)"
            )

        if self.init == "random":
            # We need to split the randomness for random initialization
            child_states = reshape(
                rng.bit_generator.spawn(n_reps * n_params * n_folds),
                (n_params, n_reps, n_folds),
            )
        else:
            # Dummy otherwise
            child_states = full((n_params, n_reps, n_folds), None)

        out = Parallel(
            n_jobs=self.n_jobs, verbose=10 if self.verbose else 0
        )(
            delayed(_penalized_cv_score)(
                structure_penalty,
                max_rank,
                factor_penalty,
                vs,
                ds,
                us,
                train_indices,
                test_indices,
                child_states[idx_param, idx_init, idx_fold],
            )
            for idx_param, (
                structure_penalty,
                max_rank,
                factor_penalty,
            ) in enumerate(parameter_grid)
            for idx_init, (vs, ds, us) in inits()
            for idx_fold, (train_indices, test_indices) in enumerate(
                cv_splits
            )
        )

        (
            scores,
            elapsed_process_times,
            est_max_rank,
            structural_zeros,
            factor_zeros,
        ) = zip(*out)

        for i in range(n_folds):
            results[f"{self.score}_fold{i}"] = [nan] * n_params

        best_runs = [-1] * n_params
        best_score = [inf] * n_params
        for idx_params, scores_params in enumerate(
            split(asarray(scores), n_params)
        ):
            for idx_init, scores_inits in enumerate(
                split(scores_params, n_reps)
            ):
                if mean(scores_inits) < best_score[idx_params]:
                    best_score[idx_params] = mean(scores_inits)
                    best_runs[idx_params] = idx_init
                    for i in range(n_folds):
                        results[f"{self.score}_fold{i}"][idx_params] = (
                            scores_inits[i]
                        )

        elapsed_process_times = split(
            asarray(elapsed_process_times), n_params
        )
        results["mean_elapsed_process_time"] = [
            mean(ts) for ts in elapsed_process_times
        ]
        results["std_elapsed_process_time"] = [
            std(ts) for ts in elapsed_process_times
        ]

        results["est_max_rank"] = [
            mean(
                est_max_rank[
                    (
                        idx_params * n_reps * n_folds
                        + best_runs[idx_params] * n_folds
                    ) : (
                        idx_params * n_reps * n_folds
                        + (best_runs[idx_params] + 1) * n_folds
                    )
                ]
            )
            for idx_params in range(n_params)
        ]
        results["structural_zeros"] = [
            mean(
                structural_zeros[
                    (
                        idx_params * n_reps * n_folds
                        + best_runs[idx_params] * n_folds
                    ) : (
                        idx_params * n_reps * n_folds
                        + (best_runs[idx_params] + 1) * n_folds
                    )
                ]
            )
            for idx_params in range(n_params)
        ]

        results["factor_zeros"] = [
            mean(
                factor_zeros[
                    (
                        idx_params * n_reps * n_folds
                        + best_runs[idx_params] * n_folds
                    ) : (
                        idx_params * n_reps * n_folds
                        + (best_runs[idx_params] + 1) * n_folds
                    )
                ]
            )
            for idx_params in range(n_params)
        ]

    # Post-processing on the full dictionary. Same for both cases
    scores = vstack(
        [results[f"{self.score}_fold{i}"] for i in range(n_folds)]
    )
    results.update(
        {
            f"mean_{self.score}": scores.mean(0),
            f"std_{self.score}": scores.std(0),
        }
    )

    self.cv_results_ = results

    if self.verbose:
        print("Re-fit final estimator")

    if self.refit.startswith("mean"):
        self.best_index_ = argmax(results[f"mean_{self.score}"])
    elif self.refit.startswith("1se"):
        # Choose the solution with maximal structure sparsity within
        # 1 standard error of the best solution
        max_index = argmax(results[f"mean_{self.score}"])

        candidates = flatnonzero(
            results[f"mean_{self.score}"]
            >= (
                results[f"mean_{self.score}"][max_index]
                - results[f"std_{self.score}"][max_index]
            )
        )

        # Primarily choose the solution with the most
        # structural zeros and then select the solution with the
        # most factor zeros if factor sparsity was requested
        structural_zeros = asarray(
            [results["structural_zeros"][i] for i in candidates]
        )
        most_sz_candidates = candidates[
            flatnonzero(structural_zeros == max(structural_zeros))
        ]

        factor_zeros = [
            results["factor_zeros"][i] for i in most_sz_candidates
        ]

        self.best_index_ = most_sz_candidates[argmax(factor_zeros)]

    structure_penalty, max_rank, factor_penalty = parameter_grid[
        self.best_index_
    ]

    if self.verbose:
        print(
            "Best fit with\n"
            f"  structure_penalty = {structure_penalty}\n"
            f"  max_rank = {max_rank}\n"
            f"  factor_penalty = {factor_penalty}\n\n"
            "  estimated max_rank = "
            f"{results['est_max_rank'][self.best_index_]}"
        )

    # Re-fit best run on all data
    if self.cv_strategy == "structure_first_debiased_cv":
        # Load respective penalized estimator from cache
        est = load(
            tmppath
            / f"{self.best_index_}_{best_runs[self.best_index_]}.pkl"
        )

        self.best_estimator_: SolrCMF = clone(base_est)

        if self.refit.endswith("debiased"):
            self.best_estimator_.set_params(
                init="custom",
                init_kwargs={"reduce_max_rank": True},
                factor_pruning=False,  # Set to False always
            )
            self.best_estimator_.fit(
                X,
                structure_pattern=est.structure_pattern(),
                factor_pattern=est.factor_pattern(),
                vs=est.vs_,
                ds=est.ds_,
                us=est.us_ if hasattr(est, "us_") else None,
            )
        elif self.refit.endswith("penalized"):
            self.best_estimator_.set_params(
                structure_penalty=structure_penalty,
                max_rank=max_rank,
                factor_penalty=factor_penalty,
                init="custom",
            )
            self.best_estimator_.fit(
                X,
                vs=est.vs_,
                ds=est.ds_,
                us=est.us_ if hasattr(est, "us_") else None,
            )

        tmpdir.cleanup()
    elif self.cv_strategy == "penalized_cv":
        # A penalized fit needs to be performed irrespectively.
        # Either because it is the final fit or because we need the
        # structure/factor pattern.
        if self.init == "custom":
            assert vs is not None and ds is not None
            vs_init = vs[best_runs[self.best_index_]]
            ds_init = ds[best_runs[self.best_index_]]
            if us is not None:
                us_init = us[best_runs[self.best_index_]]
            else:
                us_init = None
        else:
            vs_init = None
            ds_init = None
            us_init = None

        final_est: SolrCMF = clone(base_est)
        final_est.set_params(
            structure_penalty=structure_penalty,
            max_rank=max_rank,
            factor_penalty=factor_penalty,
        )
        final_est.fit(
            X,
            vs=vs_init,
            ds=ds_init,
            us=us_init,
        )

        if self.refit.endswith("debiased"):
            final_est_debiased: SolrCMF = clone(base_est)
            final_est_debiased.set_params(
                init="custom",
                init_kwargs={"reduce_max_rank": True},
                factor_pruning=False,  # Set to False always
            )
            final_est_debiased.fit(
                X,
                structure_pattern=final_est.structure_pattern(),
                factor_pattern=final_est.factor_pattern(),
                vs=final_est.vs_,
                ds=final_est.ds_,
                us=final_est.us_ if hasattr(final_est, "us_") else None,
            )

            self.best_estimator_ = final_est_debiased
        elif self.refit.endswith("penalized"):
            self.best_estimator_ = final_est

    self.best_max_rank_ = self.best_estimator_.est_max_rank_

    return self

Initialization

solrcmf.multiview_init(xs, max_rank)

Compute a decomposition using the SVD of the concatenated matrices.

Parameters:

Name	Type	Description	Default
xs	dict[ViewDesc, NDArray[float64]]	Input data	required
max_rank	int	Maximum rank of the decomposition	required

Returns:

Type	Description
dict[Entity, NDArray[float64]]	A tuple (vs, ds) containing the factor matrices in vs and the
dict[ViewDesc, NDArray[float64]]	singular values in ds.

Source code in src/solrcmf/initstrategies.py

def multiview_init(
    xs: dict[ViewDesc, NDArray[float64]],
    max_rank: int,
) -> tuple[dict[Entity, NDArray[float64]], dict[ViewDesc, NDArray[float64]]]:
    """Compute a decomposition using the SVD of the concatenated matrices.

    Args:
        xs: Input data
        max_rank: Maximum rank of the decomposition

    Returns:
        A tuple (vs, ds) containing the factor matrices in vs and the
        singular values in ds.

    """
    layout = list(xs.keys())
    if len({k[0] for k in layout}) == 1:
        x_joint = hstack([x for x in xs.values()]).T
        jx = 0
        ix = 1
    elif len({k[1] for k in layout}) == 1:
        x_joint = vstack([x for x in xs.values()])
        jx = 1
        ix = 0
    else:
        raise ValueError("'xs' does not follow a multiview layout")

    u, _, vt = svd(x_joint)

    vs = {layout[0][jx]: vt.T[:, :max_rank]}

    current = 0
    for k, x in xs.items():
        vs.update({k[ix]: u[current : current + x.shape[ix], :max_rank]})
        current += x.shape[ix]

    ds = {k: diag(vs[k[0]].T @ x @ vs[k[1]]) for k, x in xs.items()}

    return vs, ds

solrcmf.best_random_init(xs, max_rank, *, n_inits=1, n_jobs=-1, rng=None, **kwargs)

Generate best unpenalized solution from random starting points.

Parameters:

Name	Type	Description	Default
xs	dict[ViewDesc, NDArray[float64]]	Input data	required
max_rank	int	Maximum rank	required
n_inits	int	Number of random starting points to test	1
n_jobs	int	Number of jobs to run concurrently, use as in joblib.Parallel	-1
rng	Generator | int | None	Random number generator (numpy.random.Generator), random seed, or None to choose the default random number generator.	None
**kwargs		Additional arguments passed to the SolrCMF estimator.	{}

Returns:

Name	Type	Description
fit	SolrCMF	The solution found with minimal objective value among all solutions obtained from the n_inits random starting points.

Source code in src/solrcmf/initstrategies.py

def best_random_init(
    xs: dict[ViewDesc, NDArray[float64]],
    max_rank: int,
    *,
    n_inits: int = 1,
    n_jobs: int = -1,
    rng: Generator | int | None = None,
    **kwargs,
) -> SolrCMF:
    """Generate best unpenalized solution from random starting points.

    Args:
        xs: Input data
        max_rank: Maximum rank
        n_inits: Number of random starting points to test
        n_jobs: Number of jobs to run concurrently, use as in [joblib.Parallel](https://joblib.readthedocs.io/en/stable/generated/joblib.Parallel.html#joblib.Parallel)
        rng: Random number generator ([numpy.random.Generator](https://numpy.org/doc/stable/reference/random/generator.html)),
            random seed, or `None` to choose the default random number
            generator.
        **kwargs: Additional arguments passed to the SolrCMF estimator.

    Returns:
        fit (SolrCMF): The solution found with minimal objective value among
            all solutions obtained from the `n_inits` random starting points.

    """
    if n_inits <= 0:
        raise ValueError("'n_init' needs to be a positive integer")

    rng = default_rng(rng)

    def init_run(
        xs: dict[ViewDesc, NDArray[float64]], rng: Generator
    ) -> SolrCMF:
        return SolrCMF(
            structure_penalty=0.0,
            max_rank=max_rank,
            factor_pruning=False,
            init="random",
            init_kwargs={"rng": rng},
            **kwargs,
        ).fit(xs)

    rng_inits = rng.spawn(n_inits)

    ests_init: list[SolrCMF] = Parallel(n_jobs=n_jobs, return_as="list")(
        delayed(init_run)(xs, ri) for ri in rng_inits
    )

    best_obj = inf

    for i in range(n_inits):
        if ests_init[i].objective_value_ < best_obj:
            best_obj = ests_init[i].objective_value_
            best_est_init = ests_init[i]

    return best_est_init

Synthetic data generation

solrcmf.simulate

Functions to simulate synthetic data.

This module provides functions to simulate synthetic data.

Classes:

Name	Description
SimulationResult	Return type of simulate.

Functions:

Name	Description
simulate	Simulate synthetic data confirming to the SolrCMF model.

SimulationResult

Return type of simulate.

Source code in src/solrcmf/simulate.py

class SimulationResult(TypedDict):
    """Return type of simulate."""

    xs_truth: dict[ViewDesc, NDArray[float64]]
    xs: dict[ViewDesc, NDArray[float64]]
    vs: dict[Entity, NDArray[float64]]

simulate(*, viewdims, factor_scales, scales=None, snr=1.0, factor_sparsity=None, rng=None)

Simulate synthetic data confirming to the SolrCMF model.

Parameters:

Name	Type	Description	Default
viewdims	Mapping[Entity, int]	A mapping of views to view dimensions.	required
factor_scales	Mapping[ViewDesc, ArrayLike]	A mapping of view descriptors to scalars or 1D arrays describing the strength of each factor.	required
scales	Mapping[ViewDesc, float] | None	A mapping of view descriptors to positive scalars which scale the factor_scales of the corresponding view descriptor.	None
snr	Mapping[ViewDesc, float] | float	A mapping of view descriptors to signal-to-noise ratios.	1.0
factor_sparsity	Mapping[Entity, float] | None	None if factors should be simulated without sparsity and a mapping of views to sparsity proportions otherwise.	None
rng	Generator | None	A random number generator or None to use the default random number generator.	None

Returns:

Type	Description
SimulationResult	A dictionary containing the following keys "xs_truth": A dictionary of view descriptors to groundtruth (noise-less) data. "xs": A dictionary of view descriptors to the simulated data. "vs": A dictionary of views to the the groundtruth factors.

Source code in src/solrcmf/simulate.py

def simulate(
    *,
    viewdims: Mapping[Entity, int],
    factor_scales: Mapping[ViewDesc, ArrayLike],
    scales: Mapping[ViewDesc, float] | None = None,
    snr: Mapping[ViewDesc, float] | float = 1.0,
    factor_sparsity: Mapping[Entity, float] | None = None,
    rng: Generator | None = None,
) -> SimulationResult:
    """Simulate synthetic data confirming to the SolrCMF model.

    Args:
        viewdims: A mapping of views to view dimensions.
        factor_scales: A mapping of view descriptors to scalars or
            1D arrays describing the strength of each factor.
        scales: A mapping of view descriptors to positive scalars
            which scale the factor_scales of the corresponding
            view descriptor.
        snr: A mapping of view descriptors to signal-to-noise ratios.
        factor_sparsity: `None` if factors should be simulated without
            sparsity and a mapping of views to sparsity proportions
            otherwise.
        rng: A random number generator or `None` to use the default
            random number generator.

    Returns:
        A dictionary containing the following keys

            - "xs_truth": A dictionary of view descriptors to groundtruth
            (noise-less) data.
            - "xs": A dictionary of view descriptors to the simulated
            data.
            - "vs": A dictionary of views to the the groundtruth factors.

    """
    if rng is None:
        rng = default_rng()

    factor_scales_ = {k: atleast_1d(v) for k, v in factor_scales.items()}
    shapes = [s.shape for s in factor_scales_.values()]
    if not all([len(s) == 1 and s == shapes[0] for s in shapes]):
        raise ValueError(
            "Each value in 'factor_scales' needs to be of shape (max_rank,)"
        )
    max_rank = shapes[0][0]
    if not all(len(k) >= 2 for k in factor_scales_.keys()):
        raise ValueError(
            "Each key in 'factor_scales' needs to be a tuple of two"
            " or more integers"
        )

    views = set([k[i] for k in factor_scales_.keys() for i in range(2)])
    if views != viewdims.keys():
        raise ValueError(
            "The keys of 'viewdims' need to appear in the first two entries"
            " of the keys of 'factor_scales'"
        )

    if scales is None:
        scales = {k: 1.0 for k in factor_scales_.keys()}

    if scales.keys() != factor_scales_.keys():
        raise ValueError(
            "'scales' needs to be compatible with 'factor_scales'"
        )
    if not all(s > 0.0 for s in scales.values()):
        raise ValueError("Each value in 'scales' needs to be positive")

    if isinstance(snr, (int, float)):
        if snr <= 0.0:
            raise ValueError("'snr' needs to be positive")
        snr = {k: float(snr) for k in factor_scales_.keys()}
    else:
        if snr.keys() != factor_scales_.keys():
            raise ValueError(
                "'snr' needs to be compatible with 'factor_scales'"
            )
        if not all(s > 0.0 for s in snr.values()):
            raise ValueError("Each value in 'snr' needs to be positive")
        snr = dict(snr)

    if factor_sparsity is None:
        vs = {
            k: qr(rng.standard_normal((p, max_rank))).Q
            for k, p in viewdims.items()
        }
    else:
        if not (
            len(factor_sparsity) == len(views)
            and factor_sparsity.keys() == views
        ):
            raise ValueError(
                "'factor_sparsity' needs to be provided for each view"
            )

        vs = {
            k: _sparse_v(p, max_rank, factor_sparsity[k], rng)
            for k, p in viewdims.items()
        }

    xs_truth = {
        k: scales[k] * vs[k[0]] @ diag(d) @ vs[k[1]].T
        for k, d in factor_scales_.items()
    }

    xs = {
        k: x
        + sqrt(sum(x**2) / (snr[k] * x.size))
        * rng.standard_normal(size=x.shape)
        for k, x in xs_truth.items()
    }

    return {
        "xs_truth": xs_truth,
        "xs": xs,
        "vs": vs,
    }

Auxiliary methods

solrcmf.LowRankImputation

Low-rank matrix imputation via alternating ridge regression.

Factorises an incomplete matrix X as U @ V.T, where U and V are estimated by alternating ridge regression on the observed entries. Missing values (NaN) are excluded from all computations and can be recovered from the fitted factors via U_ @ V_.T.

Attributes:

Name	Type	Description
U_	NDArray[float64]	Left factor matrix of shape (n_samples, max_rank).
V_	NDArray[float64]	Right factor matrix of shape (n_features, max_rank).
converged_	bool	Whether the convergence criterion was met.
n_iter_	int	Number of iterations performed.
loss_	float	Final objective value.

Methods:

Name	Description
__init__	Initialize LowRankImputation.
fit	Fit the low-rank factorisation to X.

Source code in src/solrcmf/lrimpute.py

class LowRankImputation(BaseEstimator):
    """Low-rank matrix imputation via alternating ridge regression.

    Factorises an incomplete matrix X as U @ V.T, where U and V are
    estimated by alternating ridge regression on the observed entries.
    Missing values (NaN) are excluded from all computations and can be
    recovered from the fitted factors via U_ @ V_.T.

    Attributes:
        U_ (NDArray[float64]): Left factor matrix of shape
            (n_samples, max_rank).
        V_ (NDArray[float64]): Right factor matrix of shape
            (n_features, max_rank).
        converged_ (bool): Whether the convergence criterion was met.
        n_iter_ (int): Number of iterations performed.
        loss_ (float): Final objective value.

    """

    _parameter_constraints = {
        "penalty": [Interval(Real, 0, None, closed="left")],
        "max_rank": [Interval(Integral, 1, None, closed="left")],
        "init": [StrOptions({"random", "custom"})],
        "warm_start": ["boolean"],
        "max_iter": [Interval(Integral, 1, None, closed="left")],
        "tol": [Interval(Real, 0, None, closed="left")],
        "random_state": ["random_state"],
    }

    def __init__(
        self,
        *,
        penalty: float = 1.0,
        max_rank: int = 10,
        init: str = "random",
        warm_start: bool = False,
        max_iter: int = 1000,
        tol: float = 1e-6,
        random_state: int | RandomState | None = None,
    ):
        """Initialize LowRankImputation.

        Args:
            penalty: Ridge regularisation weight applied to both U and V.
            max_rank: Number of latent factors.
            init: Initialisation strategy. "random" draws U and V from a
                standard normal distribution; "custom" uses the U and V
                provided to fit.
            warm_start: If True, reuse U_ and V_ from a previous fit as the
                starting point instead of reinitialising.
            max_iter: Maximum number of alternating update iterations.
            tol: Convergence tolerance; stops when the relative decrease in
                loss falls below tol.
            random_state: Seed or RandomState instance for reproducible random
                initialisation.

        """
        self.penalty = penalty
        self.max_rank = max_rank
        self.init = init
        self.warm_start = warm_start
        self.max_iter = max_iter
        self.tol = tol
        self.random_state = random_state

    def _more_tags(self):
        return {"allow_nan": True}

    def fit(self, X, y=None, *, U=None, V=None) -> "LowRankImputation":
        """Fit the low-rank factorisation to X.

        Args:
            X: Input matrix of shape (n_samples, n_features), possibly
                containing NaN for missing entries.
            y: Ignored.
            U: Initial left factor matrix of shape (n_samples, max_rank).
                Used when init='custom' or warm_start=True.
            V: Initial right factor matrix of shape (n_features, max_rank).
                Used when init='custom' or warm_start=True.

        Returns:
            self

        """
        self._validate_params()

        X = check_array(
            X, dtype=[float64, float32], ensure_all_finite="allow-nan"
        )

        if U is not None or V is not None:
            if self.init != "custom":
                warn(
                    "When init!='custom', provided U or V are ignored. Set"
                    " init='custom' to use them as initialization."
                )

        if self.warm_start and hasattr(self, "U_") and hasattr(self, "V_"):
            U, V = _validate_init(X, self.U_, self.V_, self.max_rank)
        elif self.init == "random":
            U, V = _random_init(X, self.max_rank, self.random_state)
        else:  # init == "custom"
            U, V = _validate_init(X, U, V, self.max_rank)

        penalty = self.penalty
        max_iter = self.max_iter
        tol = self.tol

        loss_old = _compute_loss(X, U, V, penalty)

        converged = False
        for i in range(max_iter):
            # We will solve
            # min_{u, v} 0.5 sum_{i, j obs.} (x^(i, j) - u^(i, :) v^(j, :))^2
            #            + lambda / 2 * ||u||_F^2
            #            + lambda / 2 * ||v||_F^2

            # Given fixed v this is a ridge regression problem for each
            # u^(i, :) for a subset of the rows of v
            for r in range(X.shape[0]):
                indices = flatnonzero(1 - isnan(X[r, :]))
                A = V[indices, :].T @ V[indices, :]
                fill_diagonal(A, diagonal(A) + penalty)
                b = V[indices, :].T @ X[r, :][indices]
                U[r, :] = solve(A, b)

            # Given fixed u this is a ridge regression problem for each
            # v^(j, :) for a subset of the rows of u
            for c in range(X.shape[1]):
                indices = flatnonzero(1 - isnan(X[:, c]))
                A = U[indices, :].T @ U[indices, :]
                fill_diagonal(A, diagonal(A) + penalty)
                b = U[indices, :].T @ X[:, c][indices]
                V[c, :] = solve(A, b)

            loss = _compute_loss(X, U, V, penalty)

            if (loss_old - loss) < tol * loss_old:
                converged = True
                break

            loss_old = loss

        self.converged_ = converged
        self.U_ = U
        self.V_ = V
        self.n_iter_ = i + 1
        self.n_features_in_ = X.shape[1]
        self.loss_ = loss

        return self

__init__(*, penalty=1.0, max_rank=10, init='random', warm_start=False, max_iter=1000, tol=1e-06, random_state=None)

Initialize LowRankImputation.

Parameters:

Name	Type	Description	Default
penalty	float	Ridge regularisation weight applied to both U and V.	1.0
max_rank	int	Number of latent factors.	10
init	str	Initialisation strategy. "random" draws U and V from a standard normal distribution; "custom" uses the U and V provided to fit.	'random'
warm_start	bool	If True, reuse U_ and V_ from a previous fit as the starting point instead of reinitialising.	False
max_iter	int	Maximum number of alternating update iterations.	1000
tol	float	Convergence tolerance; stops when the relative decrease in loss falls below tol.	1e-06
random_state	int | RandomState | None	Seed or RandomState instance for reproducible random initialisation.	None

Source code in src/solrcmf/lrimpute.py

def __init__(
    self,
    *,
    penalty: float = 1.0,
    max_rank: int = 10,
    init: str = "random",
    warm_start: bool = False,
    max_iter: int = 1000,
    tol: float = 1e-6,
    random_state: int | RandomState | None = None,
):
    """Initialize LowRankImputation.

    Args:
        penalty: Ridge regularisation weight applied to both U and V.
        max_rank: Number of latent factors.
        init: Initialisation strategy. "random" draws U and V from a
            standard normal distribution; "custom" uses the U and V
            provided to fit.
        warm_start: If True, reuse U_ and V_ from a previous fit as the
            starting point instead of reinitialising.
        max_iter: Maximum number of alternating update iterations.
        tol: Convergence tolerance; stops when the relative decrease in
            loss falls below tol.
        random_state: Seed or RandomState instance for reproducible random
            initialisation.

    """
    self.penalty = penalty
    self.max_rank = max_rank
    self.init = init
    self.warm_start = warm_start
    self.max_iter = max_iter
    self.tol = tol
    self.random_state = random_state

fit(X, y=None, *, U=None, V=None)

Fit the low-rank factorisation to X.

Parameters:

Name	Type	Description	Default
X		Input matrix of shape (n_samples, n_features), possibly containing NaN for missing entries.	required
y		Ignored.	None
U		Initial left factor matrix of shape (n_samples, max_rank). Used when init='custom' or warm_start=True.	None
V		Initial right factor matrix of shape (n_features, max_rank). Used when init='custom' or warm_start=True.	None

Returns:

Type	Description
LowRankImputation	self

Source code in src/solrcmf/lrimpute.py

def fit(self, X, y=None, *, U=None, V=None) -> "LowRankImputation":
    """Fit the low-rank factorisation to X.

    Args:
        X: Input matrix of shape (n_samples, n_features), possibly
            containing NaN for missing entries.
        y: Ignored.
        U: Initial left factor matrix of shape (n_samples, max_rank).
            Used when init='custom' or warm_start=True.
        V: Initial right factor matrix of shape (n_features, max_rank).
            Used when init='custom' or warm_start=True.

    Returns:
        self

    """
    self._validate_params()

    X = check_array(
        X, dtype=[float64, float32], ensure_all_finite="allow-nan"
    )

    if U is not None or V is not None:
        if self.init != "custom":
            warn(
                "When init!='custom', provided U or V are ignored. Set"
                " init='custom' to use them as initialization."
            )

    if self.warm_start and hasattr(self, "U_") and hasattr(self, "V_"):
        U, V = _validate_init(X, self.U_, self.V_, self.max_rank)
    elif self.init == "random":
        U, V = _random_init(X, self.max_rank, self.random_state)
    else:  # init == "custom"
        U, V = _validate_init(X, U, V, self.max_rank)

    penalty = self.penalty
    max_iter = self.max_iter
    tol = self.tol

    loss_old = _compute_loss(X, U, V, penalty)

    converged = False
    for i in range(max_iter):
        # We will solve
        # min_{u, v} 0.5 sum_{i, j obs.} (x^(i, j) - u^(i, :) v^(j, :))^2
        #            + lambda / 2 * ||u||_F^2
        #            + lambda / 2 * ||v||_F^2

        # Given fixed v this is a ridge regression problem for each
        # u^(i, :) for a subset of the rows of v
        for r in range(X.shape[0]):
            indices = flatnonzero(1 - isnan(X[r, :]))
            A = V[indices, :].T @ V[indices, :]
            fill_diagonal(A, diagonal(A) + penalty)
            b = V[indices, :].T @ X[r, :][indices]
            U[r, :] = solve(A, b)

        # Given fixed u this is a ridge regression problem for each
        # v^(j, :) for a subset of the rows of u
        for c in range(X.shape[1]):
            indices = flatnonzero(1 - isnan(X[:, c]))
            A = U[indices, :].T @ U[indices, :]
            fill_diagonal(A, diagonal(A) + penalty)
            b = U[indices, :].T @ X[:, c][indices]
            V[c, :] = solve(A, b)

        loss = _compute_loss(X, U, V, penalty)

        if (loss_old - loss) < tol * loss_old:
            converged = True
            break

        loss_old = loss

    self.converged_ = converged
    self.U_ = U
    self.V_ = V
    self.n_iter_ = i + 1
    self.n_features_in_ = X.shape[1]
    self.loss_ = loss

    return self

solrcmf.bicenter(X, tol=1e-16, max_iter=10)

Bicenter the input matrix allowing for missing values.

Instead of simply centering all elements around a total mean value, this function models the data as X = m + rm + cm + Y, where m is the total mean (shape ()), rm are row means (shape (n, 1)), cm are column means (shape (1, m)), and Y are residuals (shape (n, m)).

Implements the centering algorithm described in

Hastie et al. (2015) Matrix completion and low-rank SVD via fast alternating least squares. Journal of Machine Learning Research, 16(104):3367--3402, 2015.

Parameters:

Name	Type	Description	Default
X	NDArray[floating[Any]]	The input matrix	required
tol	float	Convergence tolerance	1e-16
max_iter	int	Maximum number of iterations to perform.	10

Returns:

Type	Description
tuple[NDArray[floating[Any]], floating[Any], NDArray[floating[Any]], NDArray[floating[Any]]]	A tuple (Y, m, rm, cm) containing the bi-centered matrix Y, the overall mean m, as well as row-means rm and column-means cm.

Source code in src/solrcmf/preprocess.py

def bicenter(
    X: NDArray[floating[Any]], tol: float = 1e-16, max_iter: int = 10
) -> tuple[
    NDArray[floating[Any]],
    floating[Any],
    NDArray[floating[Any]],
    NDArray[floating[Any]],
]:
    """Bicenter the input matrix allowing for missing values.

    Instead of simply centering all elements around a total mean value, this
    function models the data as `X = m + rm + cm + Y`, where m is the total
    mean (shape ()), rm are row means (shape (n, 1)), cm are column means
    (shape (1, m)), and Y are residuals (shape (n, m)).

    Implements the centering algorithm described in
    > Hastie et al. (2015) Matrix completion and low-rank SVD via fast
    > alternating least squares. Journal of Machine Learning Research,
    > 16(104):3367--3402, 2015.

    Args:
        X: The input matrix
        tol: Convergence tolerance
        max_iter: Maximum number of iterations to perform.

    Returns:
        A tuple (Y, m, rm, cm) containing the bi-centered matrix Y, the
            overall mean m, as well as row-means rm and column-means cm.

    """
    X = check_array(X, ensure_all_finite="allow-nan")
    if tol <= 0:
        raise ValueError(f"{tol=} needs to be positive")
    if not (isinstance(max_iter, Integral) and max_iter > 0):
        raise ValueError(f"{max_iter=} needs to be a positive integer")

    n, p = X.shape

    mask = logical_not(isnan(X))
    indices = flatnonzero(mask)
    row_indices = [flatnonzero(mask[i, :]) for i in range(n)]
    col_indices = [flatnonzero(mask[:, i]) for i in range(p)]

    # Initialization
    total_mean = mean(X.flat[indices])
    row_means = array(
        [
            mean(X[i, :].flat[idx]) if len(idx) > 0 else 0.0
            for i, idx in enumerate(row_indices)
        ]
    )[:, None]
    col_means = array(
        [
            mean(X[:, i].flat[idx]) if len(idx) > 0 else 0.0
            for i, idx in enumerate(col_indices)
        ]
    )[None, :]

    # Iterate
    for it in range(max_iter):
        total_mean = mean(
            X.flat[indices] - (row_means + col_means).flat[indices]
        )
        row_means = array(
            [
                (
                    mean(
                        X[i, :].flat[idx] - (total_mean + col_means).flat[idx]
                    )
                    if len(idx) > 0
                    else 0.0
                )
                for i, idx in enumerate(row_indices)
            ]
        )[:, None]
        col_means = array(
            [
                (
                    mean(
                        X[:, i].flat[idx] - (total_mean + row_means).flat[idx]
                    )
                    if len(idx) > 0
                    else 0.0
                )
                for i, idx in enumerate(col_indices)
            ]
        )[None, :]
        r_crit = _residual(
            X,
            indices,
            row_indices,
            col_indices,
            total_mean,
            row_means,
            col_means,
        )

        if r_crit <= tol:
            break

    if it + 1 == max_iter:
        warn(f"Bi-centering did not converge in {max_iter} iterations")

    Y = X.copy()
    Y[mask] -= (total_mean + row_means + col_means)[mask]

    return Y, total_mean, row_means, col_means

solrcmf.nanscale(X, scale)

Scale all non-nan values in an array.

Parameters:

Name	Type	Description	Default
X	NDArray[floating[Any]]	Input array to be scaled, possibly containing numpy.nan values.	required
scale	float	Positive scale parameter.	required

Returns:

Type	Description
NDArray[floating[Any]]	A scaled version of the input array.

Source code in src/solrcmf/preprocess.py

def nanscale(
    X: NDArray[floating[Any]], scale: float
) -> NDArray[floating[Any]]:
    """Scale all non-nan values in an array.

    Args:
        X: Input array to be scaled, possibly containing numpy.nan values.
        scale: Positive scale parameter.

    Returns:
        A scaled version of the input array.

    """
    X = check_array(X, ensure_all_finite="allow-nan")
    if scale <= 0.0:
        raise ValueError(f"{scale=} is required to be positive")

    Y = full(X.shape, nan, dtype=X.dtype)
    divide(
        X,
        scale,
        out=Y,
        where=logical_not(isnan(X)),
    )

    return Y