markowitz.core

markowitz.core

Mathematical core for Markowitz mean-variance optimization.

This subpackage provides:

• Numerically stable linear algebra primitives (linalg).
• Closed-form Merton scalars and analytic frontier (merton_scalars, frontier).
• An immutable Portfolio value object (portfolio).
• A taxonomy of domain-specific exceptions (exceptions).

AnalyticFrontier(mu: npt.ArrayLike, Sigma: npt.ArrayLike)

Closed-form mean-variance frontier for (mu, Sigma).

Construction is the only point at which Sigma is factored; all subsequent queries (gmv, tangency, efficient_return, frontier) are O(n) once the cached basis vectors are available.

Source code in src/markowitz/core/frontier.py

def __init__(self, mu: npt.ArrayLike, Sigma: npt.ArrayLike) -> None:
    mu_arr = np.asarray(mu, dtype=np.float64)
    Sigma_arr = np.asarray(Sigma, dtype=np.float64)

    if mu_arr.ndim != 1:
        raise NumericalError(f"mu must be 1-D, got shape {mu_arr.shape!r}")
    if Sigma_arr.ndim != 2 or Sigma_arr.shape[0] != Sigma_arr.shape[1]:
        raise NumericalError(f"Sigma must be square 2-D, got shape {Sigma_arr.shape!r}")
    if mu_arr.shape[0] != Sigma_arr.shape[0]:
        raise NumericalError(
            "Dimension mismatch between mu and Sigma: "
            f"{mu_arr.shape[0]} vs {Sigma_arr.shape[0]}"
        )
    if not np.all(np.isfinite(mu_arr)):
        raise NumericalError("mu contains non-finite values")

    # Validate that Sigma is symmetric positive definite.  This
    # raises SingularCovarianceError on failure, matching the
    # documented contract.
    psd_check(Sigma_arr)

    n = mu_arr.shape[0]
    ones = np.ones(n, dtype=np.float64)
    rhs = np.column_stack((ones, mu_arr))
    sol = cholesky_solve(Sigma_arr, rhs)
    Sinv_one = sol[:, 0]
    Sinv_mu = sol[:, 1]

    A = float(ones @ Sinv_one)
    B = float(ones @ Sinv_mu)
    C = float(mu_arr @ Sinv_mu)
    D = A * C - B * B

    self._mu = mu_arr
    self._Sigma = Sigma_arr
    self._Sinv_one = Sinv_one
    self._Sinv_mu = Sinv_mu
    self._A = A
    self._B = B
    self._C = C
    self._D = D
    self._n = int(mu_arr.shape[0])

abcd: MertonABCD property

The Merton scalars associated with (mu, Sigma).

is_degenerate: bool property

True when D is essentially zero.

A degenerate frontier collapses to a single point in mean- variance space because mu is collinear with 1 (every asset has the same expected return). Efficient-return queries are still well-defined at the GMV expected return but are infeasible for any other target.

n_assets: int property

Number of assets in the investment universe.

efficient_return(mu_p: float) -> Portfolio

Return the minimum-variance portfolio with expected return mu_p.

Source code in src/markowitz/core/frontier.py

def efficient_return(self, mu_p: float) -> Portfolio:
    """Return the minimum-variance portfolio with expected return ``mu_p``."""

    from markowitz.core.portfolio import Portfolio  # noqa: PLC0415 - circular-import guard

    mu_target = float(mu_p)
    if self.is_degenerate:
        raise InfeasibleFrontierError(
            "Frontier is degenerate (D ~ 0); efficient_return is "
            "undefined for arbitrary target returns."
        )

    lam1 = (self._A * mu_target - self._B) / self._D
    lam2 = (self._C - self._B * mu_target) / self._D
    w = lam1 * self._Sinv_mu + lam2 * self._Sinv_one
    var = (self._A * mu_target * mu_target - 2.0 * self._B * mu_target + self._C) / self._D
    # Floor to zero to absorb sub-eps rounding at the GMV vertex.
    var = max(var, 0.0)
    return Portfolio(
        weights=w,
        expected_return=mu_target,
        volatility=float(np.sqrt(var)),
        sharpe=None,
    )

frontier(n_points: int = 100, *, mu_min: float | None = None, mu_max: float | None = None) -> list[Portfolio]

Return n_points portfolios spanning [mu_min, mu_max].

Defaults straddle the GMV expected return symmetrically: if no bounds are supplied, the frontier extends from mu_GMV to mu_GMV + 2 * (max(mu) - mu_GMV). Callers that need a specific range should pass it explicitly.

Source code in src/markowitz/core/frontier.py

def frontier(
    self,
    n_points: int = 100,
    *,
    mu_min: float | None = None,
    mu_max: float | None = None,
) -> list[Portfolio]:
    """Return ``n_points`` portfolios spanning ``[mu_min, mu_max]``.

    Defaults straddle the GMV expected return symmetrically: if no
    bounds are supplied, the frontier extends from ``mu_GMV`` to
    ``mu_GMV + 2 * (max(mu) - mu_GMV)``.  Callers that need a
    specific range should pass it explicitly.
    """

    if n_points < 2:
        raise ValueError(f"n_points must be >= 2, got {n_points}")

    mu_gmv = self._B / self._A
    if mu_min is None:
        mu_min = float(mu_gmv)
    if mu_max is None:
        top = float(np.max(self._mu))
        mu_max = float(mu_gmv + 2.0 * (top - mu_gmv))
        if mu_max <= mu_min:
            mu_max = mu_min + 1.0

    if mu_max < mu_min:
        raise ValueError(f"mu_max ({mu_max}) must be >= mu_min ({mu_min})")

    grid = np.linspace(float(mu_min), float(mu_max), int(n_points))
    return [self.efficient_return(float(m)) for m in grid]

gmv() -> Portfolio

Return the global minimum-variance portfolio.

Source code in src/markowitz/core/frontier.py

def gmv(self) -> Portfolio:
    """Return the global minimum-variance portfolio."""

    from markowitz.core.portfolio import Portfolio  # noqa: PLC0415 - circular-import guard

    w = self._Sinv_one / self._A
    var = 1.0 / self._A
    mu_p = self._B / self._A
    return Portfolio(
        weights=w,
        expected_return=float(mu_p),
        volatility=float(np.sqrt(var)),
        sharpe=None,
    )

tangency(rf: float) -> Portfolio

Return the tangency portfolio for risk-free rate rf.

Raises:

Type	Description
InfeasibleFrontierError	If rf coincides with the GMV expected return (i.e. B - A * rf is numerically zero), making the tangency line parallel to the frontier.

Source code in src/markowitz/core/frontier.py

def tangency(self, rf: float) -> Portfolio:
    """Return the tangency portfolio for risk-free rate ``rf``.

    Raises
    ------
    InfeasibleFrontierError
        If ``rf`` coincides with the GMV expected return (i.e.
        ``B - A * rf`` is numerically zero), making the tangency
        line parallel to the frontier.
    """

    from markowitz.core.portfolio import Portfolio  # noqa: PLC0415 - circular-import guard

    rf_f = float(rf)
    denom = self._B - self._A * rf_f
    scale = max(abs(self._B), self._A * abs(rf_f), 1.0)
    if abs(denom) < 1e-12 * scale:
        raise InfeasibleFrontierError(
            "Tangency portfolio is undefined: risk-free rate coincides "
            f"with mu_GMV = {self._B / self._A:.10g} (B - A*rf = {denom:.3e})"
        )

    # w = Sigma^-1 (mu - rf * 1) / (B - A * rf)
    # Reuse cached basis vectors instead of a second Cholesky solve.
    w = (self._Sinv_mu - rf_f * self._Sinv_one) / denom
    excess = self._C - 2.0 * self._B * rf_f + self._A * rf_f * rf_f
    # Algebraically excess >= 0 for a PD Sigma, but floor it to guard
    # against floating-point chatter near zero.
    excess = max(excess, 0.0)
    sharpe = float(np.sqrt(excess))
    er = float(w @ self._mu)
    var = float(w @ self._Sigma @ w)
    vol = float(np.sqrt(max(var, 0.0)))
    return Portfolio(
        weights=w,
        expected_return=er,
        volatility=vol,
        sharpe=sharpe,
    )

CoreError

Bases: Exception

Base class for all errors raised by markowitz.core.

InfeasibleFrontierError

Bases: CoreError

Raised when the requested frontier point does not exist.

Examples include requesting the tangency portfolio at a risk-free rate that coincides with the GMV expected return, or asking for an efficient portfolio at a return that lies outside the attainable band of a degenerate frontier.

MertonABCD

Bases: NamedTuple

The four Merton scalars characterising an analytic frontier.

NumericalError

Bases: CoreError

Raised for numerical failures that are not covered by the more specific subclasses (e.g. NaN inputs, non-finite intermediate results, or breakdowns of an iterative solver).

Portfolio(weights: FloatArray, expected_return: float, volatility: float, sharpe: float | None = None) dataclass

A snapshot of a portfolio's weights and analytic performance.

Attributes:

Name	Type	Description
weights	FloatArray	Length-n float64 array. Stored as a read-only view; the caller-provided array is copied to guarantee immutability of the :class:Portfolio instance.
expected_return	float	w^T mu evaluated at construction time.
volatility	float	sqrt(w^T Sigma w) evaluated at construction time.
sharpe	float | None	Sharpe ratio relative to whatever risk-free rate the caller chose; None if undefined.

variance() -> float

Return the portfolio variance volatility ** 2.

Source code in src/markowitz/core/portfolio.py

def variance(self) -> float:
    """Return the portfolio variance ``volatility ** 2``."""

    return self.volatility * self.volatility

SingularCovarianceError(msg: str = '', *, condition_number: float = float('inf'), min_eigenvalue: float = 0.0)

Bases: CoreError

Raised when a covariance matrix is not positive definite.

The instance attaches diagnostic information that is helpful for downstream regularization heuristics or user-facing error messages.

Parameters:

Name	Type	Description	Default
msg	str	Human-readable explanation of the failure.	''
condition_number	float	Ratio lambda_max / lambda_min of the offending matrix. Set to float('inf') when the matrix is exactly singular.	float('inf')
min_eigenvalue	float	Smallest eigenvalue of the offending matrix (possibly negative).	0.0

Source code in src/markowitz/core/exceptions.py

def __init__(
    self,
    msg: str = "",
    *,
    condition_number: float = float("inf"),
    min_eigenvalue: float = 0.0,
) -> None:
    super().__init__(msg)
    self.condition_number = condition_number
    self.min_eigenvalue = min_eigenvalue

cholesky_solve(Sigma: npt.ArrayLike, b: npt.ArrayLike, *, check_finite: bool = True, overwrite_b: bool = False) -> FloatArray

Solve Sigma x = b for a symmetric positive definite Sigma.

The Cholesky factorization is preferred over forming the explicit inverse: it is twice as fast for SPD systems and avoids the catastrophic cancellation associated with inverting near-singular matrices.

Parameters:

Name	Type	Description	Default
Sigma	ArrayLike	Symmetric positive definite matrix of shape (n, n).	required
b	ArrayLike	Right-hand side of shape (n,) or (n, k).	required
check_finite	bool	Forwarded to SciPy; when True, the inputs are screened for NaN/inf before factorization.	True
overwrite_b	bool	Permit SciPy to overwrite b to save memory. The default of False is the safe choice for general callers.	False

Returns:

Name	Type	Description
x	FloatArray	Solution array with the same shape as b.

Raises:

Type	Description
SingularCovarianceError	If Sigma is not positive definite.
NumericalError	If check_finite is true and an input contains NaN/inf.

Source code in src/markowitz/core/linalg.py

def cholesky_solve(
    Sigma: npt.ArrayLike,
    b: npt.ArrayLike,
    *,
    check_finite: bool = True,
    overwrite_b: bool = False,
) -> FloatArray:
    """Solve ``Sigma x = b`` for a symmetric positive definite ``Sigma``.

    The Cholesky factorization is preferred over forming the explicit
    inverse: it is twice as fast for SPD systems and avoids the
    catastrophic cancellation associated with inverting near-singular
    matrices.

    Parameters
    ----------
    Sigma:
        Symmetric positive definite matrix of shape ``(n, n)``.
    b:
        Right-hand side of shape ``(n,)`` or ``(n, k)``.
    check_finite:
        Forwarded to SciPy; when ``True``, the inputs are screened for
        ``NaN``/``inf`` before factorization.
    overwrite_b:
        Permit SciPy to overwrite ``b`` to save memory.  The default of
        ``False`` is the safe choice for general callers.

    Returns
    -------
    x:
        Solution array with the same shape as ``b``.

    Raises
    ------
    SingularCovarianceError
        If ``Sigma`` is not positive definite.
    NumericalError
        If ``check_finite`` is true and an input contains ``NaN``/``inf``.
    """

    Sigma_arr = np.asarray(Sigma, dtype=np.float64)
    b_arr = np.asarray(b, dtype=np.float64)

    if Sigma_arr.ndim != 2 or Sigma_arr.shape[0] != Sigma_arr.shape[1]:
        raise ValueError(f"Sigma must be a square 2-D matrix, got shape {Sigma_arr.shape!r}")
    if b_arr.shape[0] != Sigma_arr.shape[0]:
        raise ValueError(f"Shape mismatch: Sigma is {Sigma_arr.shape}, b is {b_arr.shape}")

    if check_finite:
        if not np.all(np.isfinite(Sigma_arr)):
            raise NumericalError("Sigma contains non-finite values")
        if not np.all(np.isfinite(b_arr)):
            raise NumericalError("b contains non-finite values")

    try:
        cho = cho_factor(Sigma_arr, lower=True, check_finite=False)
    except np.linalg.LinAlgError as exc:
        # SciPy raises LinAlgError when the matrix is not PD.
        eigs = eigvalsh(Sigma_arr) if Sigma_arr.size else np.array([0.0])
        raise SingularCovarianceError(
            "Cholesky factorization failed; matrix is not positive definite",
            condition_number=float("inf"),
            min_eigenvalue=float(np.min(eigs)) if eigs.size else 0.0,
        ) from exc

    solution: FloatArray = cho_solve(cho, b_arr, overwrite_b=overwrite_b, check_finite=False)
    return solution

compute_abcd(mu: npt.ArrayLike, Sigma: npt.ArrayLike) -> MertonABCD

Compute Merton's (A, B, C, D) scalars.

Parameters:

Name	Type	Description	Default
mu	ArrayLike	Expected-return vector of length n.	required
Sigma	ArrayLike	Symmetric positive definite covariance matrix of shape (n, n).	required

Returns:

Name	Type	Description
A	class:`MertonABCD` named tuple.

Raises:

Type	Description
NumericalError	If shapes are inconsistent or mu contains non-finite values.
SingularCovarianceError	Propagated from :func:cholesky_solve if Sigma is not PD.

Source code in src/markowitz/core/merton_scalars.py

def compute_abcd(mu: npt.ArrayLike, Sigma: npt.ArrayLike) -> MertonABCD:
    """Compute Merton's ``(A, B, C, D)`` scalars.

    Parameters
    ----------
    mu:
        Expected-return vector of length ``n``.
    Sigma:
        Symmetric positive definite covariance matrix of shape ``(n, n)``.

    Returns
    -------
    A :class:`MertonABCD` named tuple.

    Raises
    ------
    NumericalError
        If shapes are inconsistent or ``mu`` contains non-finite values.
    SingularCovarianceError
        Propagated from :func:`cholesky_solve` if ``Sigma`` is not PD.
    """

    mu_arr = np.asarray(mu, dtype=np.float64)
    Sigma_arr = np.asarray(Sigma, dtype=np.float64)

    if mu_arr.ndim != 1:
        raise NumericalError(f"mu must be a 1-D vector, got shape {mu_arr.shape!r}")
    if Sigma_arr.ndim != 2 or Sigma_arr.shape[0] != Sigma_arr.shape[1]:
        raise NumericalError(f"Sigma must be a square 2-D matrix, got shape {Sigma_arr.shape!r}")
    if mu_arr.shape[0] != Sigma_arr.shape[0]:
        raise NumericalError(
            f"Dimension mismatch: mu has length {mu_arr.shape[0]} but Sigma is {Sigma_arr.shape}"
        )
    if not np.all(np.isfinite(mu_arr)):
        raise NumericalError("mu contains non-finite values")

    n = mu_arr.shape[0]
    ones = np.ones(n, dtype=np.float64)

    # Stack the two right-hand sides so we do a single Cholesky solve
    # rather than two.  This is purely a performance optimization; the
    # numerical result is identical to two separate solves.
    rhs = np.column_stack((ones, mu_arr))
    sol = cholesky_solve(Sigma_arr, rhs)
    Sinv_one = sol[:, 0]
    Sinv_mu = sol[:, 1]

    A = float(ones @ Sinv_one)
    B = float(ones @ Sinv_mu)
    C = float(mu_arr @ Sinv_mu)
    D = A * C - B * B

    return MertonABCD(A=A, B=B, C=C, D=D)

condition_number(Sigma: npt.ArrayLike) -> float

Return lambda_max / lambda_min for Sigma.

Returns float('inf') if the matrix is not positive definite, contains non-finite entries, or has a non-positive smallest eigenvalue. No exception is raised.

Source code in src/markowitz/core/linalg.py

def condition_number(Sigma: npt.ArrayLike) -> float:
    """Return ``lambda_max / lambda_min`` for ``Sigma``.

    Returns ``float('inf')`` if the matrix is not positive definite,
    contains non-finite entries, or has a non-positive smallest
    eigenvalue.  No exception is raised.
    """

    arr = np.asarray(Sigma, dtype=np.float64)
    if arr.ndim != 2 or arr.shape[0] != arr.shape[1] or arr.size == 0:
        return float("inf")
    if not np.all(np.isfinite(arr)):
        return float("inf")
    sym = 0.5 * (arr + arr.T)
    try:
        eigs = eigvalsh(sym)
    except np.linalg.LinAlgError:  # pragma: no cover - eigvalsh rarely fails
        return float("inf")
    lam_min = float(np.min(eigs))
    lam_max = float(np.max(eigs))
    if lam_min <= 0.0:
        return float("inf")
    return lam_max / lam_min

portfolio_performance(weights: npt.ArrayLike, mu: npt.ArrayLike, Sigma: npt.ArrayLike, *, rf: float | None = None) -> Portfolio

Compute analytic performance for a given weight vector.

Parameters:

Name	Type	Description	Default
weights	ArrayLike	Length-n weight vector. No sum-to-one or sign constraints are imposed; the function simply evaluates the quadratic form.	required
mu	ArrayLike	Expected returns.	required
Sigma	ArrayLike	Covariance matrix.	required
rf	float | None	Optional risk-free rate. When supplied, sharpe is populated; otherwise it is None.	None

Returns:

Name	Type	Description
A	class:`Portfolio` instance.

Raises:

Type	Description
NumericalError	If shapes are inconsistent or inputs contain non-finite values.

Source code in src/markowitz/core/portfolio.py

def portfolio_performance(
    weights: npt.ArrayLike,
    mu: npt.ArrayLike,
    Sigma: npt.ArrayLike,
    *,
    rf: float | None = None,
) -> Portfolio:
    """Compute analytic performance for a given weight vector.

    Parameters
    ----------
    weights:
        Length-``n`` weight vector.  No sum-to-one or sign constraints
        are imposed; the function simply evaluates the quadratic form.
    mu:
        Expected returns.
    Sigma:
        Covariance matrix.
    rf:
        Optional risk-free rate.  When supplied, ``sharpe`` is
        populated; otherwise it is ``None``.

    Returns
    -------
    A :class:`Portfolio` instance.

    Raises
    ------
    NumericalError
        If shapes are inconsistent or inputs contain non-finite values.
    """

    w = np.asarray(weights, dtype=np.float64)
    mu_arr = np.asarray(mu, dtype=np.float64)
    S = np.asarray(Sigma, dtype=np.float64)

    if w.ndim != 1:
        raise NumericalError(f"weights must be 1-D, got shape {w.shape!r}")
    if mu_arr.ndim != 1 or mu_arr.shape != w.shape:
        raise NumericalError(f"mu must be 1-D with shape {w.shape}, got {mu_arr.shape!r}")
    if S.ndim != 2 or S.shape != (w.shape[0], w.shape[0]):
        raise NumericalError(f"Sigma must be {w.shape[0]}x{w.shape[0]}, got {S.shape!r}")
    if not (np.all(np.isfinite(w)) and np.all(np.isfinite(mu_arr)) and np.all(np.isfinite(S))):
        raise NumericalError("Inputs contain non-finite values")

    er = float(w @ mu_arr)
    var = float(w @ S @ w)
    # Variance can be slightly negative for numerically zero portfolios.
    # Clip to zero before sqrt to avoid spurious NaNs.
    vol = float(np.sqrt(max(var, 0.0)))

    sharpe: float | None = None if rf is None else ((er - float(rf)) / vol if vol > 0.0 else None)

    return Portfolio(weights=w, expected_return=er, volatility=vol, sharpe=sharpe)

psd_check(Sigma: npt.ArrayLike, *, tol: float = 1e-10, require_symmetric: bool = True, symmetry_tol: float = 1e-10) -> None

Validate that Sigma is symmetric positive definite.

Parameters:

Name	Type	Description	Default
Sigma	ArrayLike	Candidate covariance matrix.	required
tol	float	Lower bound for the smallest eigenvalue. Matrices whose minimum eigenvalue is below tol are rejected.	1e-10
require_symmetric	bool	Reject matrices that are not symmetric within symmetry_tol.	True
symmetry_tol	float	Absolute tolerance for the ||S - S.T||_inf check.	1e-10

Raises:

Type	Description
SingularCovarianceError	If the matrix is not symmetric (when required), not positive definite, or contains non-finite entries.
ValueError	If the matrix is not 2-D and square.

Source code in src/markowitz/core/linalg.py

def psd_check(
    Sigma: npt.ArrayLike,
    *,
    tol: float = 1e-10,
    require_symmetric: bool = True,
    symmetry_tol: float = 1e-10,
) -> None:
    """Validate that ``Sigma`` is symmetric positive definite.

    Parameters
    ----------
    Sigma:
        Candidate covariance matrix.
    tol:
        Lower bound for the smallest eigenvalue.  Matrices whose
        minimum eigenvalue is below ``tol`` are rejected.
    require_symmetric:
        Reject matrices that are not symmetric within ``symmetry_tol``.
    symmetry_tol:
        Absolute tolerance for the ``||S - S.T||_inf`` check.

    Raises
    ------
    SingularCovarianceError
        If the matrix is not symmetric (when required), not positive
        definite, or contains non-finite entries.
    ValueError
        If the matrix is not 2-D and square.
    """

    arr = np.asarray(Sigma, dtype=np.float64)
    if arr.ndim != 2 or arr.shape[0] != arr.shape[1]:
        raise ValueError(f"Sigma must be a square 2-D matrix, got shape {arr.shape!r}")
    if not np.all(np.isfinite(arr)):
        raise SingularCovarianceError(
            "Sigma contains non-finite values",
            condition_number=float("inf"),
            min_eigenvalue=float("nan"),
        )

    if require_symmetric:
        asymmetry = float(np.max(np.abs(arr - arr.T))) if arr.size else 0.0
        if asymmetry > symmetry_tol:
            raise SingularCovarianceError(
                f"Sigma is not symmetric (max |S - S.T| = {asymmetry:.3e})",
                condition_number=float("inf"),
                min_eigenvalue=float("nan"),
            )

    # Symmetric eigenvalues; cheaper and more accurate than ``eig``.
    sym = 0.5 * (arr + arr.T)
    eigs = eigvalsh(sym) if sym.size else np.array([0.0])
    lam_min = float(np.min(eigs)) if eigs.size else 0.0
    lam_max = float(np.max(eigs)) if eigs.size else 0.0
    if lam_min < tol:
        cond = float("inf") if lam_min <= 0 else lam_max / lam_min
        raise SingularCovarianceError(
            (f"Sigma is not positive definite (min eigenvalue = {lam_min:.3e} < tol = {tol:.1e})"),
            condition_number=cond,
            min_eigenvalue=lam_min,
        )

regularize(Sigma: npt.ArrayLike, *, eps: float | Literal['auto'] = 'auto', method: Literal['ridge'] = 'ridge') -> FloatArray

Return a regularized copy of Sigma that is safely SPD.

The default ridge method adds eps * I to a symmetrized copy of Sigma. When eps='auto', a scale-aware value is chosen based on the trace of Sigma and machine precision.

Parameters:

Name	Type	Description	Default
Sigma	ArrayLike	Input matrix.	required
eps	float | Literal['auto']	Either an explicit non-negative scalar or the string 'auto'.	'auto'
method	Literal['ridge']	Currently only 'ridge' is supported.	'ridge'

Returns:

Type	Description
A new ``(n, n)`` ``float64`` array.

Source code in src/markowitz/core/linalg.py

def regularize(
    Sigma: npt.ArrayLike,
    *,
    eps: float | Literal["auto"] = "auto",
    method: Literal["ridge"] = "ridge",
) -> FloatArray:
    """Return a regularized copy of ``Sigma`` that is safely SPD.

    The default ``ridge`` method adds ``eps * I`` to a symmetrized copy
    of ``Sigma``.  When ``eps='auto'``, a scale-aware value is chosen
    based on the trace of ``Sigma`` and machine precision.

    Parameters
    ----------
    Sigma:
        Input matrix.
    eps:
        Either an explicit non-negative scalar or the string ``'auto'``.
    method:
        Currently only ``'ridge'`` is supported.

    Returns
    -------
    A new ``(n, n)`` ``float64`` array.
    """

    arr = np.asarray(Sigma, dtype=np.float64)
    if arr.ndim != 2 or arr.shape[0] != arr.shape[1]:
        raise ValueError(f"Sigma must be a square 2-D matrix, got shape {arr.shape!r}")
    if method != "ridge":
        raise ValueError(f"Unsupported regularization method: {method!r}")

    n = arr.shape[0]
    sym = 0.5 * (arr + arr.T)

    if eps == "auto":
        # Scale by mean diagonal so eps is dimensionally consistent with
        # Sigma; floor at 1e-8 so that the regularized matrix is well
        # above the default ``psd_check`` tolerance (1e-10) even when
        # the input has zero or tiny entries.
        diag_mean = float(np.mean(np.abs(np.diag(sym)))) if n else 0.0
        eps_val = max(diag_mean * 1e-8, 1e-8)
    else:
        eps_val = float(eps)
        if eps_val < 0.0:
            raise ValueError(f"eps must be non-negative, got {eps_val}")

    result: np.ndarray = sym + eps_val * np.eye(n, dtype=sym.dtype)
    return result