mavats

113@dataclass
114class EnvelopeMARResult(MARResult):
115    """Fixed-dimension EMAR estimate and local-optimization diagnostics.
116
117    The model paper and references of :func:`fit_envelope_mar` also apply to
118    all inherited forecast/covariance methods. ``row_envelope`` and
119    ``column_envelope`` are orthonormal bases; compare their projectors, not
120    their coordinates. Every lag's output lies in these spaces and each
121    innovation covariance is reduced by the corresponding space.
122
123    ``optimization_history`` records every attempted sweep, with row/column
124    selected-start and per-start objectives, gradient norms and stop reasons.
125    ``converged`` requires both outer likelihood and selected inner gradient
126    criteria, not global optimality or stationarity of the time series.
127    ``warmup_converged`` concerns only the unrestricted initialization.
128    ``objective_history`` contains accepted complete EMAR sweeps on scaled
129    data, excluding the unrestricted warmup. Covariance scaling and physical
130    likelihood conversion follow MARResult. Dense diagnostic allocations are
131    guarded by ``max_dense_elements``; matrix-form forecasts avoid them.
132    """
133
134    envelope_dims: tuple = (0, 0)
135    row_envelope: np.ndarray = field(default_factory=lambda: np.empty((0, 0)))
136    column_envelope: np.ndarray = field(default_factory=lambda: np.empty((0, 0)))
137    optimization_history: tuple = ()
138    warmup_converged: bool | None = None
139    stop_reason: str = "max_iter"
140    likelihood_converged: bool = False
141    inner_converged: bool = False
142    max_dense_elements: int = 10_000_000
143
144    def _guard_dense(self, count):
145        if count > self.max_dense_elements:
146            raise ValueError(
147                "dense diagnostic exceeds max_dense_elements; matrix-form forecasting is still available"
148            )
149
150    @property
151    def coefficients(self):
152        self._guard_dense(self.order * self.intercept.size**2)
153        return super().coefficients
154
155    @property
156    def spectral_radius(self):
157        if self.order > 1:
158            self._guard_dense((self.order * self.intercept.size) ** 2)
159        return super().spectral_radius
160
161    @property
162    def log_likelihood(self):
163        """Physical-unit conditional Gaussian likelihood, without squaring scale."""
164        return float(
165            -self.objective_history[-1] - self.residuals.size * np.log(self.data_scale)
166        )
167
168    def residual_covariance(self, *, ddof=0):
169        self._guard_dense(self.intercept.size**2)
170        return super().residual_covariance(ddof=ddof)

173def fit_envelope_mar(
174    X,
175    envelope_dims,
176    *,
177    order=1,
178    fit_intercept=True,
179    max_iter=100,
180    tol=1e-8,
181    inner_max_iter=200,
182    inner_tol=1e-6,
183    inner_starts=3,
184    random_state=0,
185    initial=None,
186    initial_covariance=None,
187    initial_envelopes=None,
188    warmup_max_iter=100,
189    max_dense_elements=10_000_000,
190):
191    """Fit Gaussian envelope MAR(p) with fixed positive envelope dimensions.
192
193    Parameters
194    ----------
195    X : array_like, shape (T, m, n)
196        Finite real time-first observations. All lags share the same N=T-p
197        regression rows. Missing observations are not supported.
198    envelope_dims : (u, v)
199        Response envelope dimensions, 1 <= u <= m and 1 <= v <= n.
200        A_l = R eta_l and B_l = C xi_l; past inputs are NOT projected onto
201        R and C. Covariances have material and orthogonal-complement blocks.
202        Full dimensions give the unrestricted separable MAR likelihood.
203    order : int, default 1
204        Number of distinct bilinear lag terms with shared output envelopes.
205    fit_intercept : bool, default True
206        Profile a free intercept by centering response and each lag over the
207        same regression sample; restore the intercept in original units.
208    max_iter, tol : int, float
209        Complete row/column sweeps and relative negative-loglikelihood
210        tolerance on data divided by its maximum absolute entry.
211    inner_max_iter, inner_tol, inner_starts : int, float, int
212        Grassmann iterations, absolute gradient tolerance, and fresh starts
213        per block. Starts are the smallest residual eigenvectors, smallest
214        response eigenvectors, then random orthonormal bases. The previous
215        envelope is always an additional start. No global optimum guarantee.
216    random_state : int or numpy.random.Generator, default 0
217        Local random starts; does not change NumPy's global RNG.
218    initial : (left, right), optional
219        Finite factors shaped (p,m,m), (p,n,n); 2-D accepted for p=1.
220        Bypasses unrestricted MAR-MLE warmup. Without this, the paper's
221        unrestricted warmup uses up to warmup_max_iter sweeps.
222    initial_covariance : (row, column), optional
223        SPD covariance factors in physical units, divided by data_scale
224        individually for computation. With no initial coefficients these
225        initialize the unrestricted warmup. Otherwise default to identities
226        in scaled-data units.
227    initial_envelopes : (R, C), optional
228        Orthonormal bases shaped (m,u), (n,v), used as extra first starts.
229        Without supplied bases, an unrestricted warmup also supplies leading
230        left singular vectors of its concatenated row/column coefficients.
231        These are additional numerical starts, not fixed envelope estimates.
232    warmup_max_iter : int, default 100
233        Unrestricted MAR-MLE initialization budget, separately diagnosed.
234    max_dense_elements : int, default 10000000
235        Guard for dense result diagnostics and automatic projection warmup.
236        Larger warmups use matrix-sized identity initialization instead.
237
238    Returns
239    -------
240    EnvelopeMARResult
241        Bilinear factors, reducing spaces, separable covariance, forecasts
242        and inner/outer optimization diagnostics. Inspect ``converged``.
243
244    Notes
245    -----
246    Implements equations (15)--(16) using joint-lag, whitened SVD regressions
247    and the actual logdet Grassmann profile, not PCA or reduced-rank ALS.
248    Singular weighted designs and numerically singular residual moments
249    raise: no ridge, covariance flooring or pseudodeterminant is substituted.
250    Each nonzero lag pair is scale/sign normalized separately. Stability is
251    diagnosed, never imposed. Only fixed-dimension point estimation is
252    implemented: no dimension/order selection, sparse SEMAR, standard errors
253    or the separate one-sided identity-coefficient model.
254
255    References
256    ----------
257    Samadi, S. Y. and De Alwis, T. P. (2026), "Envelope Matrix Autoregressive
258    Models", Journal of Business & Economic Statistics 44(2), 397--412,
259    equations (13), (15)--(16), Algorithm 1,
260    https://doi.org/10.1080/07350015.2025.2537404.
261    Numerical solver and source-equation corrections are documented in
262    docs/envelope-notes.md; this is not an exact replication of paper code.
263    """
264    order = positive_int(order, "order")
265    original = as_series(X, min_samples=order + 2)
266    m, n = original.shape[1:]
267    dims = ranks_tuple(envelope_dims, (m, n))
268    if dims is None:
269        raise ValueError("envelope_dims must specify two positive dimensions")
270    max_iter = positive_int(max_iter, "max_iter")
271    inner_max_iter = positive_int(inner_max_iter, "inner_max_iter")
272    inner_starts = positive_int(inner_starts, "inner_starts")
273    warmup_max_iter = positive_int(warmup_max_iter, "warmup_max_iter")
274    max_dense_elements = positive_int(max_dense_elements, "max_dense_elements")
275    tol = finite_scalar(tol, "tol", minimum=0)
276    inner_tol = finite_scalar(inner_tol, "inner_tol", minimum=0)
277    if not isinstance(fit_intercept, (bool, np.bool_)):
278        raise ValueError("fit_intercept must be boolean")
279    rng = random_generator(random_state)
280    scale = float(np.max(np.abs(original))) or 1.0
281    data = original / scale
282    raw_lags = [data[order - lag : len(data) - lag] for lag in range(1, order + 1)]
283    mean = data[order:].mean(axis=0) if fit_intercept else np.zeros((m, n))
284    lag_means = [
285        z.mean(axis=0) if fit_intercept else np.zeros((m, n)) for z in raw_lags
286    ]
287    response = data[order:] - mean
288    lags = [z - mu for z, mu in zip(raw_lags, lag_means)]
289    bases = [None, None]
290    if initial_envelopes is not None:
291        if len(initial_envelopes) != 2:
292            raise ValueError("initial_envelopes must contain row and column bases")
293        bases = [
294            _basis(b, d, r, "initial envelope")
295            for b, d, r in zip(initial_envelopes, (m, n), dims)
296        ]
297    covariances = [np.eye(m), np.eye(n)]
298    if initial_covariance is not None:
299        if len(initial_covariance) != 2:
300            raise ValueError("initial_covariance must contain two matrices")
301        for i, (value, d) in enumerate(zip(initial_covariance, (m, n))):
302            raw = np.asarray(value)
303            if np.iscomplexobj(raw) or raw.dtype.kind not in "biuf":
304                raise ValueError("initial covariance must be real numeric")
305            cov = np.array(raw, dtype=float, copy=True) / scale
306            if cov.shape != (d, d) or not np.allclose(cov, cov.T, rtol=1e-10, atol=0):
307                raise ValueError(
308                    "initial covariance must be symmetric and dimensionally compatible"
309                )
310            covariances[i] = _spd(cov, "initial covariance")[0]
311    warmup_converged = None
312    if initial is None:
313        warmup = fit_mar(
314            data,
315            method="mle",
316            order=order,
317            fit_intercept=fit_intercept,
318            max_iter=warmup_max_iter,
319            tol=tol,
320            covariance_floor=0,
321            initial_covariance=tuple(covariances),
322            init="auto" if order * (m * n) ** 2 <= max_dense_elements else "identity",
323        )
324        left, right = warmup.left.copy(), warmup.right.copy()
325        covariances = [warmup.row_covariance.copy(), warmup.column_covariance.copy()]
326        warmup_converged = warmup.converged
327        if initial_envelopes is None:
328            # Covariance-eigenvector starts alone can be trapped in immaterial
329            # directions when material variance is high. Use the training-only
330            # warmup's coefficient spaces as an additional first start.
331            bases = [
332                np.linalg.svd(np.concatenate(coef, axis=1), full_matrices=False)[0][
333                    :, :rank
334                ]
335                for coef, rank in zip((left, right), dims)
336            ]
337    else:
338        if len(initial) != 2:
339            raise ValueError("initial must contain left and right coefficients")
340        coefficients = []
341        for value, d in zip(initial, (m, n)):
342            raw = np.asarray(value)
343            if np.iscomplexobj(raw) or raw.dtype.kind not in "biuf":
344                raise ValueError("initial coefficients must be real numeric")
345            array = np.array(raw, dtype=float, copy=True)
346            if order == 1 and array.ndim == 2:
347                array = array[None]
348            if array.shape != (order, d, d) or not np.isfinite(array).all():
349                raise ValueError(
350                    "initial coefficient dimensions must match order and X"
351                )
352            coefficients.append(array)
353        left, right = coefficients
354    row, column = covariances
355    history, optimization = [], []
356    likelihood_converged = inner_converged = False
357    reason = "max_iter"
358    options = dict(starts=inner_starts, rng=rng, max_iter=inner_max_iter, tol=inner_tol)
359    transposed_lags = [z.transpose(0, 2, 1) for z in lags]
360    for iteration in range(1, max_iter + 1):
361        a, new_row, r, rd = _block_update(
362            response, lags, right, column, dims[0], initial=bases[0], **options
363        )
364        new_left = np.stack(np.split(a, order, axis=1))
365        b, new_column, c, cd = _block_update(
366            response.transpose(0, 2, 1),
367            transposed_lags,
368            new_left,
369            new_row,
370            dims[1],
371            initial=bases[1],
372            **options,
373        )
374        new_right = np.stack(np.split(b, order, axis=1))
375        for lag in range(order):
376            # Never replace a zero A with an out-of-envelope coordinate vector.
377            if np.any(new_left[lag]):
378                new_left[lag], new_right[lag] = _normalize(
379                    new_left[lag], new_right[lag]
380                )
381        norm = np.linalg.norm(new_row)
382        new_row, new_column = new_row / norm, new_column * norm
383        residual = response - _prediction(lags, new_left, new_right)
384        value = _nll(residual, new_row, new_column)
385        if not np.isfinite(value):
386            raise FloatingPointError("envelope likelihood exceeds floating-point range")
387        accepted = not history or value <= history[-1] + 1e-10 * max(
388            1, abs(history[-1])
389        )
390        optimization.append(
391            dict(
392                iteration=iteration,
393                row=rd,
394                column=cd,
395                accepted=accepted,
396                objective=value,
397            )
398        )
399        if not accepted:
400            reason = "likelihood_increase"
401            likelihood_converged = False
402            break
403        likelihood_converged = bool(
404            history and abs(history[-1] - value) <= tol * max(1, abs(history[-1]))
405        )
406        inner_converged = rd["converged"] and cd["converged"]
407        history.append(value)
408        left, right, row, column, bases = (
409            new_left,
410            new_right,
411            new_row,
412            new_column,
413            [r, c],
414        )
415        if likelihood_converged and inner_converged:
416            reason = "likelihood_and_gradient_tolerance"
417            break
418    intercept = mean.copy()
419    for a, mu, b in zip(left, lag_means, right):
420        intercept -= a @ mu @ b.T
421    fitted = mean + _prediction(lags, left, right)
422    with np.errstate(over="ignore", invalid="ignore"):
423        physical = [
424            intercept * scale,
425            fitted * scale,
426            (data[order:] - fitted) * scale,
427            row * scale,
428            column * scale,
429        ]
430    if not all(np.isfinite(value).all() for value in physical):
431        raise FloatingPointError(
432            "envelope fit exceeds floating-point range in physical units"
433        )
434    return EnvelopeMARResult(
435        left=left,
436        right=right,
437        intercept=physical[0],
438        fitted_values=physical[1],
439        residuals=physical[2],
440        method="envelope-mle",
441        converged=likelihood_converged and inner_converged,
442        n_iter=len(history),
443        objective_history=np.asarray(history),
444        row_covariance=physical[3],
445        column_covariance=physical[4],
446        data_scale=scale,
447        _history=original[-order:].copy(),
448        envelope_dims=dims,
449        row_envelope=bases[0],
450        column_envelope=bases[1],
451        optimization_history=tuple(optimization),
452        warmup_converged=warmup_converged,
453        stop_reason=reason,
454        likelihood_converged=likelihood_converged,
455        inner_converged=inner_converged,
456        max_dense_elements=max_dense_elements,
457    )

 94@dataclass
 95class RankSelectionStep:
 96    """One mode's criterion evaluated at its sequential update.
 97
 98    ``eigenvalues * exp(log_eigenvalue_scale)`` are the physical fourth-power
 99    spectrum; neither exponentiation nor physical spectra need be representable.
100    ``log_penalty`` is in those same physical units. IC ``scores`` share the
101    positive factor ``exp(log_score_scale)``; ER scores are ordinary ratios.
102    ``log_scores`` retain physical IC values or dimensionless ER log ratios.
103    Rounding may tie displayed scores even when stable criterion comparisons
104    distinguish them. The exact-tie convention is the smallest candidate rank.
105    ``projected_shape`` records the data dimensions used for this moment, NOT
106    the original dimensions retained in the theoretical penalty.
107
108    References
109    ----------
110    Han, Chen and Zhang (2022), Rank Determination in Tensor Factor Model,
111    equations (1)--(2), (7)--(8). https://doi.org/10.1214/22-EJS1991
112    """
113
114    rank: int
115    eigenvalues: np.ndarray
116    log_eigenvalue_scale: float
117    log_penalty: float
118    candidate_ranks: np.ndarray
119    scores: np.ndarray
120    log_scores: np.ndarray
121    log_score_scale: float
122    projected_shape: tuple

133@dataclass
134class TensorRankSelectionResult:
135    """Selected dynamic ranks, fitted Tucker spaces, and complete sweep history.
136
137    ``history[0]`` is the unprojected selection; subsequent entries contain one
138    sequential sweep. ``initial_ranks`` are those unprojected choices, while
139    ``starting_ranks`` are the conservative projection initialization. The final
140    ``ranks`` are per-mode criterion choices, not necessarily minimal multilinear
141    ranks: any zero mode implies an identically zero centered signal, regardless
142    of other modes' reported ranks. ``signal_ranks`` makes this distinction.
143    Convergence means ranks AND loading projectors stopped changing, not correct
144    rank recovery. Iterative spectra refer to the corresponding sequential
145    update, not a new simultaneous recomputation using all final loadings.
146
147    References
148    ----------
149    Han, Chen and Zhang (2022), Rank Determination in Tensor Factor Model,
150    Section 3.2 and Remarks 5--7. https://doi.org/10.1214/22-EJS1991
151    """
152
153    ranks: tuple
154    loadings: tuple
155    factors: np.ndarray
156    signal: np.ndarray
157    residuals: np.ndarray
158    mean: np.ndarray
159    history: tuple
160    initial_ranks: tuple
161    starting_ranks: tuple
162    rank_history: tuple
163    projector_changes: tuple
164    method: str
165    criterion: str
166    penalty: int
167    nu: float
168    c0: float
169    log_penalty_multipliers: tuple
170    max_ranks: tuple
171    lags: tuple
172    data_scale: float
173    n_iter: int
174    converged: bool
175    stop_reason: str
176
177    @property
178    def diagnostics(self):
179        return self.history[-1]
180
181    @property
182    def eigenvalues(self):
183        return tuple(step.eigenvalues for step in self.diagnostics)
184
185    @property
186    def signal_ranks(self):
187        return (0,) * len(self.ranks) if 0 in self.ranks else self.ranks
188
189    def transform(self, X):
190        """Project new data using only fitted spaces and the training mean.
191
192        References
193        ----------
194        Han, Chen and Zhang (2022), Rank Determination in Tensor Factor Model,
195        equations (4)--(5). https://doi.org/10.1214/22-EJS1991
196        Fixed-loading Tucker projection is reconstruction, not forecasting.
197        """
198        X = as_series(X, min_samples=1, ndim=len(self.ranks) + 1)
199        if X.shape[1:] != self.mean.shape:
200            raise ValueError("observation dimensions must match the fitted model")
201        axes = tuple(range(1, X.ndim))
202        scale = np.maximum(
203            np.max(np.abs(X), axis=axes, keepdims=True),
204            float(np.max(np.abs(self.mean))),
205        )
206        scale = np.where(scale > 0, scale, 1.0)
207        return _restore(
208            _project(X / scale - self.mean / scale, self.loadings),
209            scale,
210            "transformed factors",
211        )
212
213    def inverse_transform(self, factors):
214        """Reconstruct with the training mean, including empty zero-rank cores.
215
216        References
217        ----------
218        Han, Chen and Zhang (2022), Rank Determination in Tensor Factor Model,
219        equation (4). https://doi.org/10.1214/22-EJS1991
220        Restoring a fitted deterministic mean is the centering convention.
221        """
222        raw = np.asarray(factors)
223        if (
224            np.iscomplexobj(raw)
225            or raw.dtype.kind not in "biuf"
226            or raw.ndim != len(self.ranks) + 1
227            or raw.shape[0] < 1
228            or raw.shape[1:] != self.ranks
229        ):
230            raise ValueError("factors must be finite real cores with the fitted ranks")
231        factors = np.asarray(raw, dtype=float)
232        if not np.isfinite(factors).all():
233            raise ValueError("factors must be finite real cores with the fitted ranks")
234        axes = tuple(range(1, factors.ndim))
235        scale = np.maximum(
236            np.max(np.abs(factors), axis=axes, keepdims=True, initial=0),
237            float(np.max(np.abs(self.mean))),
238        )
239        scale = np.where(scale > 0, scale, 1.0)
240        return _restore(
241            _project(factors / scale, self.loadings, inverse=True) + self.mean / scale,
242            scale,
243            "reconstructed signal",
244        )

266def select_tensor_rank(
267    X,
268    *,
269    method="tipup",
270    criterion="er",
271    penalty=1,
272    max_ranks=None,
273    lags=1,
274    iterative=False,
275    nu=0.0,
276    penalty_multiplier=1.0,
277    log_penalty_multiplier=None,
278    c0=0.1,
279    initial_ranks=None,
280    max_iter=100,
281    tol=1e-8,
282    center=True,
283):
284    """Select Tucker ranks by the published TOPUP/TIPUP IC or ER criteria.
285
286    Parameters
287    ----------
288    X : array_like, shape (T, d1, ..., dK)
289        Finite real training observations, K >= 1 and every d_k >= 2.
290    method : {'topup', 'tipup'}
291        Past-first lag outer/inner product moments, Section 3.2.
292    criterion : {'ic', 'er'}
293        Equation (1) permits 0,...,m*; equation (2) permits 1,...,m*.
294    penalty : int, 1,...,5
295        Exact g1,...,g5 (IC) or h1,...,h5 (ER), equations (7)--(8).
296    max_ranks : sequence of int, optional
297        Search bounds m*_k < d_k; defaults to floor(d_k/2). Zero bounds are
298        allowed only for IC. Full-dimension search extensions are not provided.
299    lags : int, default 1
300        h0: accumulate lags 1,...,h0, without averaging the moment over lags.
301    iterative : bool, default False
302        Sequential equation (6) updates, reselecting every rank in each sweep.
303    nu : float in [0, 1], default 0
304        IC strength tuning, ideally weakest-factor exponent delta1. Zero is the
305        paper's default strong-factor setting, not an estimate of strength.
306    penalty_multiplier, log_penalty_multiplier : scalar or sequence, optional
307        Positive multiplicative constants, one per mode if desired. The log
308        alternative supports constants outside float range and requires the
309        ordinary multiplier to remain its scalar default 1. Constants apply
310        to FIXED physical fourth-power penalties, not relative eigenvalue
311        heuristics. Under X -> a X, add 4 log(abs(a)) to log multipliers to
312        preserve the criterion. Keeping constants fixed need NOT preserve ranks.
313    c0 : positive float, default 0.1
314        Constant in h1; ignored by other displayed penalty families.
315    initial_ranks : sequence of int, optional
316        Iterative projection initialization only, up to d_k. Default is
317        min(2*r_initial, r_initial+3, d_k), following Remark 5 with dimension cap.
318    max_iter : positive int, default 100
319        Maximum sequential sweeps, excluding the unprojected initialization.
320    tol : positive float, default 1e-8
321        Stop when ranks agree and the maximum Frobenius projector change < tol.
322    center : bool, default True
323        Remove and retain the training mean, as recommended in Remark 6 when
324        demeaning does not remove a factor dimension (constant factors violate
325        that premise). Held-out transformations use this training mean only.
326
327    Notes
328    -----
329    Original d and d_k enter penalties even after projection. Only moments use
330    projected dimensions. No penalty calibration, confidence intervals, ER
331    mock-zero eigenvalue, or full-dimension rank extension is implied. Serially
332    white noise, informative lag moments, factor strengths and the rate
333    conditions of Section 4 are substantive assumptions; numerical convergence
334    does not check them. TIPUP can lose factors by lag-inner-product cancellation.
335    IC zero-mode selections have exact empty-core reconstruction. With an empty
336    projected mode subsequent moments are zero, an absorbing loss of information
337    which the positive-rank asymptotic theory does not remedy.
338
339    References
340    ----------
341    Han, Chen and Zhang (2022), Rank Determination in Tensor Factor Model,
342    equations (1), (2), (6)--(8), Remarks 5--7.
343    https://doi.org/10.1214/22-EJS1991
344    https://arxiv.org/html/2011.07131v3
345    """
346    method, criterion = str(method).lower(), str(criterion).lower()
347    if method not in {"topup", "tipup"}:
348        raise ValueError("method must be 'topup' or 'tipup'")
349    if criterion not in {"ic", "er"}:
350        raise ValueError("criterion must be 'ic' or 'er'")
351    penalty = positive_int(penalty, "penalty")
352    if penalty > 5:
353        raise ValueError("penalty must be in 1,...,5")
354    if not isinstance(iterative, (bool, np.bool_)):
355        raise ValueError("iterative must be boolean")
356    h0 = positive_int(lags, "lags")
357    max_iter = positive_int(max_iter, "max_iter")
358    tol = finite_scalar(tol, "tol", minimum=0)
359    nu = finite_scalar(nu, "nu", minimum=0)
360    c0 = finite_scalar(c0, "c0", minimum=0)
361    if tol == 0 or c0 == 0 or nu > 1:
362        raise ValueError("tol and c0 must be positive, and nu must be <= 1")
363    ndim = np.ndim(X)
364    if ndim < 2:
365        raise ValueError("X needs time and at least one spatial mode")
366    X, work, mean, scale = _prepare(X, center, ndim=ndim)
367    shape, n = X.shape[1:], len(X)
368    if min(shape) < 2 or h0 >= n:
369        raise ValueError("all mode dimensions must be >= 2 and lags < sample size")
370    if max_ranks is None:
371        max_ranks = tuple(d // 2 for d in shape)
372    else:
373        try:
374            max_ranks = tuple(max_ranks)
375        except TypeError as exc:
376            raise ValueError("max_ranks must contain one bound per mode") from exc
377        if len(max_ranks) != len(shape):
378            raise ValueError("max_ranks must contain one bound per mode")
379        max_ranks = tuple(
380            positive_int(r, "max_rank", minimum=0 if criterion == "ic" else 1)
381            for r in max_ranks
382        )
383        if any(r >= d for r, d in zip(max_ranks, shape)):
384            raise ValueError("max_ranks must be strictly below mode dimensions")
385    if log_penalty_multiplier is None:
386        log_multipliers = _per_mode(
387            penalty_multiplier, len(shape), "penalty_multiplier"
388        )
389    else:
390        if (
391            np.ndim(penalty_multiplier) != 0
392            or finite_scalar(penalty_multiplier, "penalty_multiplier") != 1
393        ):
394            raise ValueError(
395                "supply either penalty_multiplier or log_penalty_multiplier"
396            )
397        log_multipliers = _per_mode(
398            log_penalty_multiplier,
399            len(shape),
400            "log_penalty_multiplier",
401            logarithmic=True,
402        )
403    log_penalties = tuple(
404        _log_penalty(shape, n, h0, k, criterion, penalty, nu, c0) + c
405        for k, c in enumerate(log_multipliers)
406    )
407    lag_values = tuple(range(1, h0 + 1))
408    log_scale = 4 * math.log(scale)
409
410    def update(projected, mode):
411        moment = (
412            np.zeros((shape[mode], shape[mode]))
413            if projected.size == 0
414            else _lagged_moment(_unfold_series(projected, mode), lag_values, method)
415        )
416        if not np.isfinite(moment).all():
417            raise ValueError("lag moment cannot be represented after normalization")
418        values, vectors = np.linalg.eigh((moment + moment.T) * 0.5)
419        values = np.maximum(values[::-1], 0)
420        vectors = _canonical_signs(vectors[:, ::-1])
421        chosen = _criterion(
422            values, log_scale, log_penalties[mode], criterion, max_ranks[mode]
423        )
424        step = RankSelectionStep(
425            eigenvalues=values,
426            log_eigenvalue_scale=log_scale,
427            log_penalty=log_penalties[mode],
428            projected_shape=projected.shape[1:],
429            **chosen,
430        )
431        return vectors, step
432
433    first = [update(work, k) for k in range(len(shape))]
434    initial = tuple(pair[1].rank for pair in first)
435    if initial_ranks is not None:
436        if not iterative:
437            raise ValueError("initial_ranks apply only to iterative selection")
438        try:
439            starting = tuple(initial_ranks)
440        except TypeError as exc:
441            raise ValueError("initial_ranks must contain one rank per mode") from exc
442        if len(starting) != len(shape):
443            raise ValueError("initial_ranks must contain one rank per mode")
444        starting = tuple(
445            positive_int(r, "initial_rank", minimum=0 if criterion == "ic" else 1)
446            for r in starting
447        )
448        if any(r > d for r, d in zip(starting, shape)):
449            raise ValueError("initial_ranks cannot exceed mode dimensions")
450    else:
451        starting = (
452            tuple(min(2 * r, r + 3, d) for r, d in zip(initial, shape))
453            if iterative
454            else initial
455        )
456    loadings = [pair[0][:, :r] for pair, r in zip(first, starting)]
457    history = [tuple(pair[1] for pair in first)]
458    rank_history, changes = [starting], []
459    converged, iteration = not iterative, 0
460    if iterative:
461        for iteration in range(1, max_iter + 1):
462            old = list(loadings)
463            steps = []
464            for mode in range(len(shape)):
465                vectors, step = update(_project(work, loadings, skip=mode), mode)
466                loadings[mode] = vectors[:, : step.rank]
467                steps.append(step)
468            current = tuple(u.shape[1] for u in loadings)
469            change = max(
470                float(np.linalg.norm(u @ u.T - v @ v.T)) for u, v in zip(loadings, old)
471            )
472            history.append(tuple(steps))
473            changes.append(change)
474            rank_history.append(current)
475            if current == tuple(u.shape[1] for u in old) and change < tol:
476                converged = True
477                break
478    ranks = tuple(u.shape[1] for u in loadings)
479    core = _project(work, loadings)
480    reconstructed = _project(core, loadings, inverse=True) + mean / scale
481    signal = _restore(reconstructed, scale, "fitted signal")
482    factors = _restore(core, scale, "fitted factors")
483    residuals = _restore(X / scale - reconstructed, scale, "fitted residuals")
484    return TensorRankSelectionResult(
485        ranks,
486        tuple(loadings),
487        factors,
488        signal,
489        residuals,
490        mean,
491        tuple(history),
492        initial,
493        starting,
494        tuple(rank_history),
495        tuple(changes),
496        ("i" if iterative else "") + method,
497        criterion,
498        penalty,
499        nu,
500        c0,
501        log_multipliers,
502        max_ranks,
503        lag_values,
504        scale,
505        iteration,
506        converged,
507        (
508            "rank_and_projector_tolerance"
509            if iterative and converged
510            else ("max_iter" if iterative else "noniterative")
511        ),
512    )