API Reference¶

def cholesky(cov: np.ndarray, rhs: np.ndarray | None = None) -> np.ndarray:
    """Compute the upper triangular Cholesky factor, or solve a linear system.

    When called without *rhs*, returns the upper triangular factor R such that
    R.T @ R = cov.  When *rhs* is given, solves ``cov @ x = rhs`` using the
    Cholesky decomposition for numerical stability, falling back to LU
    decomposition if *cov* is not positive-definite.

    Args:
        cov: A positive definite covariance matrix of shape (n, n).
        rhs: Optional right-hand side vector or matrix. When provided the
            system ``cov @ x = rhs`` is solved and *x* is returned.

    Returns:
        The upper triangular Cholesky factor R when *rhs* is ``None``, or the
        solution x to ``cov @ x = rhs`` otherwise.

    Raises:
        np.linalg.LinAlgError: When *rhs* is ``None`` and *cov* is not
            positive-definite, or when *rhs* is given and both Cholesky and
            LU-based solves fail.

    Example:
        Decomposition (no rhs):

        >>> import numpy as np
        >>> from cvx.linalg import cholesky
        >>> cov = np.array([[4.0, 2.0], [2.0, 5.0]])
        >>> R = cholesky(cov)
        >>> np.allclose(R.T @ R, cov)
        True

        Solve (with rhs):

        >>> cholesky(np.eye(2), np.array([1.0, 2.0])).tolist()
        [1.0, 2.0]
        >>> cholesky(np.array([[4.0, 0.0], [0.0, 9.0]]), np.array([8.0, 27.0])).tolist()
        [2.0, 3.0]
    """
    if rhs is None:
        return _cholesky(cov).transpose()
    try:
        upper = _cholesky(cov).transpose()
        return np.linalg.solve(upper, np.linalg.solve(upper.T, rhs))
    except np.linalg.LinAlgError:
        return np.linalg.solve(cov, rhs)

def is_positive_definite(matrix: np.ndarray) -> bool:
    """Return True if *matrix* is symmetric positive-definite, False otherwise.

    The check is performed via an attempted Cholesky decomposition — the most
    numerically reliable way to test positive-definiteness for symmetric matrices.

    This function is side-effect-free: it raises no exceptions and emits no
    warnings.  It is suitable for use as a guard before passing a matrix to a
    linear solver.

    Args:
        matrix: Square matrix to test.

    Returns:
        ``True`` if the matrix is positive-definite, ``False`` otherwise.

    Example:
        >>> import numpy as np
        >>> from cvx.linalg import is_positive_definite
        >>> is_positive_definite(np.eye(3))
        True
        >>> is_positive_definite(np.array([[1.0, 2.0], [2.0, 1.0]]))
        False
        >>> is_positive_definite(np.array([[1.0, 0.5], [0.5, 1.0]]))
        True
    """
    try:
        cholesky(matrix)
    except np.linalg.LinAlgError:
        return False
    else:
        return True

def eigh(matrix: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Compute eigendecomposition of a real symmetric or Hermitian matrix.

    Rows and columns with non-finite diagonal entries are excluded before
    decomposition.

    Args:
        matrix: Square symmetric/Hermitian input matrix.

    Returns:
        A tuple ``(eigenvalues, eigenvectors)`` as returned by ``np.linalg.eigh``
        on the valid submatrix. Eigenvalues are sorted in ascending order.
    """
    _, submatrix = valid(matrix)
    return np.linalg.eigh(submatrix)

def eigvalsh(matrix: np.ndarray) -> np.ndarray:
    """Return eigenvalues of a real symmetric or Hermitian matrix.

    Rows and columns with non-finite diagonal entries are excluded before
    decomposition.

    Args:
        matrix: Square symmetric/Hermitian input matrix.

    Returns:
        Eigenvalues of the valid submatrix in ascending order.
    """
    eigenvalues, _ = eigh(matrix)
    return eigenvalues

def eigvals(matrix: np.ndarray) -> np.ndarray:
    """Return the eigenvalues of a square matrix.

    This routine supports general (non-symmetric) square matrices and may
    return complex eigenvalues.

    Args:
        matrix: Square input matrix.

    Returns:
        Eigenvalues of ``matrix`` as returned by ``numpy.linalg.eigvals``.

    Raises:
        NotAMatrixError: If *matrix* is not two-dimensional.
        NonSquareMatrixError: If *matrix* is not square.
    """
    if matrix.ndim != 2:
        raise NotAMatrixError(matrix.ndim)

    if matrix.shape[0] != matrix.shape[1]:
        raise NonSquareMatrixError(matrix.shape[0], matrix.shape[1])

    return np.linalg.eigvals(matrix)

def qr(matrix: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Compute the reduced QR decomposition of a 2-D matrix.

    Args:
        matrix: Input matrix with shape ``(m, n)``.

    Returns:
        A tuple ``(Q, R)`` matching ``np.linalg.qr(matrix, mode="reduced")``.

    Raises:
        DimensionMismatchError: If ``matrix`` is not two-dimensional.

    Example:
        >>> import numpy as np
        >>> from cvx.linalg import qr
        >>> q, r = qr(np.array([[1.0, 2.0], [3.0, 4.0]]))
        >>> np.allclose(q @ r, np.array([[1.0, 2.0], [3.0, 4.0]]))
        True
    """
    if matrix.ndim != 2:
        raise DimensionMismatchError(matrix.ndim, 2)

    return np.linalg.qr(matrix, mode="reduced")

def svd(matrix: Matrix) -> tuple[Matrix, Matrix, Matrix]:
    """Compute the compact singular value decomposition of a matrix.

    This is a thin wrapper around ``numpy.linalg.svd`` with
    ``full_matrices=False``.

    Args:
        matrix: Input matrix of shape ``(m, n)``.

    Returns:
        Tuple ``(u, s, vt)`` such that ``matrix == u @ np.diag(s) @ vt``.
    """
    return np.linalg.svd(matrix, full_matrices=False)

def pca(returns: Matrix, n_components: int = 10) -> PCA:
    """Compute the first n principal components for a return matrix using SVD.

    Args:
        returns: Array of asset returns with shape (n_samples, n_assets).
        n_components: Number of principal components to extract. Defaults to 10.

    Returns:
        PCA named tuple containing:
            - explained_variance: Ratio of variance explained by each component
            - factors: Factor returns (scores)
            - exposure: Factor exposures (loadings)
            - cov: Factor covariance matrix
            - systematic: Returns explained by factors
            - idiosyncratic: Residual returns

    Example:
        >>> import numpy as np
        >>> from cvx.linalg import pca
        >>> np.random.seed(42)
        >>> returns = np.random.randn(100, 10)
        >>> result = pca(returns, n_components=3)
        >>> bool(result.explained_variance[0] > result.explained_variance[1])
        True
        >>> factor_corr = np.corrcoef(result.factors.T)
        >>> bool(np.allclose(factor_corr, np.eye(3), atol=0.1))
        True
        >>> VtV = result.exposure @ result.exposure.T
        >>> bool(np.allclose(VtV, np.eye(3), atol=1e-10))
        True
        >>> all(result.explained_variance[i] >= result.explained_variance[i+1]
        ...     for i in range(len(result.explained_variance)-1))
        True
        >>> reconstructed = result.factors @ result.exposure
        >>> centered_systematic = result.systematic - returns.mean(axis=0)
        >>> bool(np.allclose(reconstructed, centered_systematic, atol=1e-10))
        True

    """
    x_mean = returns.mean(axis=0)
    x_centered = returns - x_mean

    u, s_full, vt = svd(x_centered)

    u = u[:, :n_components]
    s = s_full[:n_components]
    vt = vt[:n_components, :]

    factors: Matrix = u * s
    exposure: Matrix = vt
    explained_variance: Matrix = (s**2) / np.sum(s_full**2)
    cov: Matrix = np.cov(factors.T)
    systematic: Matrix = factors @ vt + x_mean
    idiosyncratic: Matrix = x_centered - factors @ vt

    return PCA(
        explained_variance=explained_variance,
        factors=factors,
        exposure=exposure,
        cov=cov,
        systematic=systematic,
        idiosyncratic=idiosyncratic,
    )

def solve(
    matrix: np.ndarray,
    rhs: np.ndarray,
    cond_threshold: float = _DEFAULT_COND_THRESHOLD,
) -> np.ndarray:
    """Solve a linear system restricted to the valid submatrix.

    Rows and columns with non-finite diagonal entries are excluded from the
    solve; the corresponding positions in the result are set to NaN.  Cholesky
    decomposition is attempted first for numerical stability and falls back to
    LU decomposition for non-positive-definite matrices.  When the condition
    number of the valid sub-matrix exceeds *cond_threshold*, an
    ``IllConditionedMatrixWarning`` is emitted.

    Args:
        matrix: Square coefficient matrix.
        rhs: Right-hand side vector.
        cond_threshold: Condition-number threshold above which a warning is
            emitted. Defaults to ``1e12``.

    Returns:
        A solution vector with the same shape as ``rhs``. Entries mapped to
        invalid rows or columns are returned as ``NaN``.

    Raises:
        NonSquareMatrixError: If the matrix is not square.
        DimensionMismatchError: If ``rhs`` size does not match the matrix dimension.
        SingularMatrixError: If the valid sub-matrix is singular.

    Example:
        >>> import numpy as np
        >>> from cvx.linalg import solve
        >>> solve(np.eye(2), np.array([1.0, 2.0])).tolist()
        [1.0, 2.0]

        NaN-masked entries are skipped:

        >>> matrix = np.array([[4.0, 0.0], [0.0, np.nan]])
        >>> solve(matrix, np.array([8.0, 1.0])).tolist()
        [2.0, nan]
    """
    if matrix.shape[0] != matrix.shape[1]:
        raise NonSquareMatrixError(matrix.shape[0], matrix.shape[1])

    if rhs.size != matrix.shape[0]:
        raise DimensionMismatchError(rhs.size, matrix.shape[0])

    solution = np.nan * np.ones(rhs.size)
    mask, submatrix = valid(matrix)

    if mask.any():
        _check_and_warn_condition(submatrix, cond_threshold)
        try:
            solution[mask] = _cholesky(submatrix, rhs[mask])
        except np.linalg.LinAlgError as exc:
            raise SingularMatrixError(str(exc)) from exc

    return solution

def lstsq(
    matrix: np.ndarray,
    rhs: np.ndarray,
    cond_threshold: float = _DEFAULT_COND_THRESHOLD,
) -> tuple[np.ndarray, np.ndarray, int, np.ndarray]:
    """Solve an overdetermined or underdetermined system in the least-squares sense.

    Rows where any entry in *matrix* or the corresponding entry in *rhs* is
    non-finite are excluded before solving.  The returned solution vector
    always has length equal to the number of columns in *matrix*.  When the
    effective condition number of the valid sub-matrix exceeds
    *cond_threshold*, an ``IllConditionedMatrixWarning`` is emitted.

    Args:
        matrix: Coefficient matrix of shape ``(m, n)``.
        rhs: Right-hand side vector of length ``m``.
        cond_threshold: Condition-number threshold above which a warning is
            emitted. Defaults to ``1e12``.

    Returns:
        A four-tuple ``(x, residuals, rank, sv)`` matching the convention of
        :func:`numpy.linalg.lstsq`:

        - ``x`` — least-squares solution of shape ``(n,)``.
        - ``residuals`` — sum of squared residuals; empty when the solution is
          not unique or all rows are invalid.
        - ``rank`` — effective rank of the valid sub-matrix.
        - ``sv`` — singular values of the valid sub-matrix in descending order.

    Raises:
        DimensionMismatchError: If ``rhs`` length does not match the number of
            rows in *matrix*.

    Example:
        >>> import numpy as np
        >>> from cvx.linalg import lstsq
        >>> A = np.array([[1.0, 1.0], [1.0, 2.0], [1.0, 3.0]])
        >>> b = np.array([6.0, 5.0, 7.0])
        >>> x, res, rank, sv = lstsq(A, b)
        >>> int(rank)
        2

        NaN rows are silently dropped:

        >>> A_nan = np.array([[1.0, 1.0], [np.nan, 2.0], [1.0, 3.0]])
        >>> b_nan = np.array([6.0, 5.0, 7.0])
        >>> x2, _, rank2, _ = lstsq(A_nan, b_nan)
        >>> int(rank2)
        2
    """
    if rhs.shape[0] != matrix.shape[0]:
        raise DimensionMismatchError(rhs.shape[0], matrix.shape[0])

    n_cols = matrix.shape[1]

    # Filter rows that contain any non-finite value in matrix or rhs.
    row_mask = np.isfinite(matrix).all(axis=1) & np.isfinite(rhs)
    sub_matrix = matrix[row_mask]
    sub_rhs = rhs[row_mask]

    if sub_matrix.shape[0] == 0:
        return np.full(n_cols, np.nan), np.array([]), 0, np.array([])

    x, residuals, rank, sv = np.linalg.lstsq(sub_matrix, sub_rhs, rcond=None)

    # Compute condition number from singular values.
    if sv.size > 0 and sv[-1] > 0:
        cond = float(sv[0] / sv[-1])
    elif sv.size > 0:
        cond = float("inf")
    else:
        cond = 1.0

    if cond > cond_threshold:
        warnings.warn(
            f"Matrix condition number {cond:.3e} exceeds threshold {cond_threshold:.3e}; "
            "results may be numerically unreliable.",
            IllConditionedMatrixWarning,
            stacklevel=2,
        )

    return x, residuals, rank, sv

def inv(
    matrix: np.ndarray,
    cond_threshold: float = _DEFAULT_COND_THRESHOLD,
) -> np.ndarray:
    """Invert a matrix restricted to the valid submatrix.

    Rows and columns with non-finite diagonal entries are excluded from the
    inversion; the corresponding rows and columns in the result are set to NaN.
    When the condition number of the valid sub-matrix exceeds *cond_threshold*,
    an ``IllConditionedMatrixWarning`` is emitted.

    Args:
        matrix: Square matrix to invert.
        cond_threshold: Condition-number threshold above which a warning is
            emitted. Defaults to ``1e12``.

    Returns:
        An inverted matrix with the same shape as *matrix*. Rows and columns
        mapped to invalid entries are returned as ``NaN``.

    Raises:
        NonSquareMatrixError: If the matrix is not square.
        SingularMatrixError: If the valid sub-matrix is singular.

    Example:
        >>> import numpy as np
        >>> from cvx.linalg import inv
        >>> np.allclose(inv(np.eye(2)), np.eye(2))
        True

        NaN-masked entries are skipped:

        >>> matrix = np.array([[4.0, 0.0], [0.0, np.nan]])
        >>> result = inv(matrix)
        >>> float(result[0, 0])
        0.25
        >>> bool(np.isnan(result[0, 1]) and np.isnan(result[1, 0]) and np.isnan(result[1, 1]))
        True
    """
    if matrix.shape[0] != matrix.shape[1]:
        raise NonSquareMatrixError(matrix.shape[0], matrix.shape[1])

    n = matrix.shape[0]
    result = np.full((n, n), np.nan)
    mask, submatrix = valid(matrix)

    if mask.any():
        _check_and_warn_condition(submatrix, cond_threshold)
        try:
            sub_inv = np.linalg.inv(submatrix)
        except np.linalg.LinAlgError as exc:
            raise SingularMatrixError(str(exc)) from exc

        idx = np.where(mask)[0]
        result[np.ix_(idx, idx)] = sub_inv

    return result

def norm(
    x: np.ndarray,
    ord: int | float | Literal["fro", "nuc"] | None = None,
) -> float:
    """Compute the norm of a vector or matrix, ignoring non-finite entries.

    Non-finite entries (NaN, inf) are treated as zero before computing the norm,
    so they contribute nothing to the result. Supports all ``ord`` values accepted
    by ``np.linalg.norm``.

    Args:
        x: Input array (1-D vector or 2-D matrix).
        ord: Order of the norm. See ``np.linalg.norm`` for valid values.
            Common choices: ``None`` (default 2-norm for vectors, Frobenius for
            matrices), ``1``, ``2``, ``np.inf``, ``'fro'``, ``'nuc'``.

    Returns:
        The norm as a float.

    Example:
        >>> import numpy as np
        >>> from cvx.linalg import norm
        >>> norm(np.array([3.0, np.nan, 4.0]))
        5.0
        >>> norm(np.array([[1.0, np.nan], [np.nan, 1.0]]), ord='fro')
        1.4142135623730951
    """
    return float(np.linalg.norm(np.where(np.isfinite(x), x, 0.0), ord=ord))

def a_norm(vector: np.ndarray, matrix: np.ndarray | None = None) -> float:
    """Calculate the generalized norm of a vector with respect to a matrix.

    Args:
        vector: The input vector.
        matrix: Optional square matrix defining the quadratic form.

    Returns:
        The Euclidean norm of the finite vector entries, or ``sqrt(v.T @ A @ v)``
        after dropping rows and columns whose diagonal entries are not finite.

    Raises:
        NonSquareMatrixError: If the matrix is not square.
        DimensionMismatchError: If the vector length does not match the matrix dimension.

    Example:
        >>> import numpy as np
        >>> from cvx.linalg import a_norm
        >>> a_norm(np.array([3.0, 4.0]))
        5.0
    """
    if matrix is None:
        return float(np.linalg.norm(vector[np.isfinite(vector)], 2))

    if matrix.shape[0] != matrix.shape[1]:
        raise NonSquareMatrixError(matrix.shape[0], matrix.shape[1])

    if vector.size != matrix.shape[0]:
        raise DimensionMismatchError(vector.size, matrix.shape[0])

    mask, submatrix = valid(matrix)
    if mask.any():
        filtered_vector = vector[mask]
        return float(np.sqrt(filtered_vector @ submatrix @ filtered_vector))

    return float("nan")

def inv_a_norm(
    vector: np.ndarray,
    matrix: np.ndarray | None = None,
    cond_threshold: float = _DEFAULT_COND_THRESHOLD,
) -> float:
    """Calculate the inverse A-norm of a vector using an optional matrix.

    If ``matrix`` is ``None``, returns the Euclidean norm of finite entries.
    Otherwise computes ``sqrt(v.T @ A^{-1} @ v)`` on the valid submatrix,
    attempting Cholesky decomposition first for numerical stability and
    falling back to LU decomposition for non-positive-definite matrices.
    When the condition number of the valid sub-matrix exceeds *cond_threshold*,
    an ``IllConditionedMatrixWarning`` is emitted.

    Args:
        vector: The input vector.
        matrix: Optional square matrix defining the quadratic form.
        cond_threshold: Condition-number threshold above which a warning is
            emitted. Defaults to ``1e12``.

    Returns:
        The Euclidean norm of the finite vector entries, or
        ``sqrt(v.T @ A^{-1} @ v)`` after dropping rows and columns whose
        diagonal entries are not finite. Returns ``nan`` when no valid entries
        exist.

    Raises:
        NonSquareMatrixError: If the matrix is not square.
        DimensionMismatchError: If the vector length does not match the matrix dimension.
        SingularMatrixError: If the valid sub-matrix is singular.

    Example:
        >>> import numpy as np
        >>> from cvx.linalg import inv_a_norm
        >>> inv_a_norm(np.array([3.0, 4.0]))
        5.0
    """
    if matrix is None:
        return float(np.linalg.norm(vector[np.isfinite(vector)], 2))

    if matrix.shape[0] != matrix.shape[1]:
        raise NonSquareMatrixError(matrix.shape[0], matrix.shape[1])

    if vector.size != matrix.shape[0]:
        raise DimensionMismatchError(vector.size, matrix.shape[0])

    mask, submatrix = valid(matrix)
    if mask.any():
        _check_and_warn_condition(submatrix, cond_threshold)
        filtered_vector = vector[mask]
        try:
            solved = _cholesky(submatrix, filtered_vector)
        except np.linalg.LinAlgError as exc:
            raise SingularMatrixError(str(exc)) from exc
        return float(np.sqrt(np.dot(filtered_vector, solved)))

    return float("nan")

def cond(matrix: np.ndarray, p: int | float | Literal["fro", "nuc"] | None = None) -> float:
    """Return the condition number of a matrix.

    Returns ``nan`` if the matrix contains any non-finite (NaN or inf) entries.
    Otherwise delegates to :func:`numpy.linalg.cond`.

    Args:
        matrix: Input matrix.
        p: Order of the norm used to compute the condition number.
            Accepts the same values as :func:`numpy.linalg.cond`
            (``None``, ``1``, ``-1``, ``2``, ``-2``, ``numpy.inf``,
            ``-numpy.inf``, ``'fro'``).  Defaults to ``None`` which
            corresponds to the 2-norm (largest singular value divided by
            the smallest).

    Returns:
        The condition number as a ``float``, or ``nan`` when the matrix
        contains non-finite entries.

    Examples:
        >>> import numpy as np
        >>> cond(np.eye(3))
        1.0
        >>> import math
        >>> math.isnan(cond(np.array([[float('nan'), 1.0], [1.0, 2.0]])))
        True
        >>> cond(np.diag([1.0, 1e10]), p=1)
        10000000000.0
    """
    if not np.all(np.isfinite(matrix)):
        return float("nan")
    return float(np.linalg.cond(matrix, p=p))

def det(
    matrix: np.ndarray,
    cond_threshold: float = _DEFAULT_COND_THRESHOLD,
) -> float:
    """Return the determinant of a square matrix.

    Rows and columns with non-finite diagonal entries are excluded before the
    computation; when no valid rows or columns remain the function returns
    ``nan``.  When the condition number of the valid sub-matrix exceeds
    *cond_threshold*, an ``IllConditionedMatrixWarning`` is emitted.

    Args:
        matrix: Square input matrix.
        cond_threshold: Condition-number threshold above which a warning is
            emitted. Defaults to ``1e12``.

    Returns:
        The determinant of the valid sub-matrix, or ``nan`` when no valid
        entries exist.

    Raises:
        NonSquareMatrixError: If the matrix is not square.

    Example:
        >>> import numpy as np
        >>> from cvx.linalg import det
        >>> det(np.eye(3))
        1.0

        NaN-masked entries are skipped:

        >>> matrix = np.array([[2.0, 0.0], [0.0, np.nan]])
        >>> det(matrix)
        2.0
    """
    if matrix.shape[0] != matrix.shape[1]:
        raise NonSquareMatrixError(matrix.shape[0], matrix.shape[1])

    mask, submatrix = valid(matrix)

    if not mask.any():
        return _SENTINEL

    _check_and_warn_condition(submatrix, cond_threshold)
    return float(np.linalg.det(submatrix))

def ewm_covariance(
    data: pl.DataFrame,
    assets: list[str],
    index_col: str,
    window: int = 30,
    is_halflife: bool = False,
    warmup: int = 0,
) -> dict[Hashable, np.ndarray]:
    """Compute the exponentially weighted covariance matrix of returns.

    EWM covariance uses the identity
    ``Cov(X, Y) = EWM(X*Y) - EWM(X)*EWM(Y)`` applied to the
    *common non-null observations* of each pair, which is equivalent
    to ``pandas.DataFrame.ewm(span).cov(bias=True)``.

    Each date is included in the result as long as at least one
    matrix entry is non-NaN.  Cells involving a late-starting asset
    are ``NaN`` until that asset has enough observations; the date is
    never dropped on account of a single asset being unavailable.
    Dates where every cell is NaN (before the warmup period is met
    for any asset) are omitted.

    Args:
        data: Polars DataFrame containing the index column and asset columns.
        assets: Ordered list of asset column names.
        index_col: Name of the index (e.g. date) column in *data*.
        window: Span (default) or half-life (when *is_halflife* is
            ``True``) of the exponential decay.  Defaults to ``30``.
        is_halflife: When ``True`` *window* is interpreted as the
            half-life; otherwise it is the EWMA span.  Defaults to
            ``False``.
        warmup: Minimum number of common observations required before
            a pair's cell is non-NaN.  Defaults to ``0`` (cells are
            non-NaN from the first shared observation).

    Returns:
        Dictionary keyed by index value (date or integer) mapping to
        a square symmetric ``numpy.ndarray`` of shape ``(n, n)``
        where ``n`` is the number of assets.  Row/column order
        matches *assets*.  Unavailable cells are ``NaN``.

    """
    if isinstance(warmup, bool) or not isinstance(warmup, int):
        raise TypeError
    if warmup < 0:
        raise NegativeWarmupError

    n = len(assets)
    min_samples = 1 if warmup == 0 else warmup

    def _ewm(expr: pl.Expr) -> pl.Expr:
        """Apply EWM mean with the configured span or half-life."""
        if is_halflife:
            return expr.ewm_mean(half_life=window, min_samples=min_samples)
        return expr.ewm_mean(span=window, min_samples=min_samples)

    cov_exprs = [
        (
            _ewm(pl.col(a) * pl.col(b))
            - _ewm(pl.when(pl.col(b).is_null()).then(None).otherwise(pl.col(a)))
            * _ewm(pl.when(pl.col(a).is_null()).then(None).otherwise(pl.col(b)))
        ).alias(f"{a}_{b}")
        for i, a in enumerate(assets)
        for b in assets[i:]
    ]

    pair_df = data.with_columns(cov_exprs).drop(assets)
    all_keys = pair_df[index_col].to_list()
    pair_arr = pair_df.drop(index_col).to_numpy()

    ii, jj = np.triu_indices(n)
    cube = np.full((len(all_keys), n, n), np.nan)
    cube[:, ii, jj] = pair_arr
    cube[:, jj, ii] = pair_arr

    has_data = ~np.all(np.isnan(cube), axis=(1, 2))
    return {k: cube[t] for t, k in enumerate(all_keys) if has_data[t]}

def rand_cov(n: int, seed: int | None = None) -> np.ndarray:
    """Construct a random positive semi-definite covariance matrix of size n x n.

    The matrix is constructed as A^T @ A where A is a random n x n matrix with
    elements drawn from a standard normal distribution. This ensures the result
    is symmetric and positive semi-definite.

    Args:
        n: Size of the covariance matrix (n x n).
        seed: Random seed for reproducibility. If None, uses the current
            random state.

    Returns:
        A random positive semi-definite n x n covariance matrix.

    Example:
        Generate a reproducible random covariance matrix:

        >>> import numpy as np
        >>> from cvx.linalg import rand_cov
        >>> cov1 = rand_cov(3, seed=42)
        >>> cov2 = rand_cov(3, seed=42)
        >>> np.allclose(cov1, cov2)
        True

        Verify positive definiteness via Cholesky decomposition:

        >>> cov = rand_cov(5, seed=123)
        >>> # If Cholesky succeeds without error, matrix is positive definite
        >>> L = np.linalg.cholesky(cov)
        >>> bool(np.allclose(L @ L.T, cov))
        True

        Eigenvalue verification:

        >>> cov = rand_cov(3, seed=99)
        >>> eigenvalues = np.linalg.eigvalsh(cov)
        >>> # All eigenvalues should be positive for PD matrix
        >>> bool(np.all(eigenvalues > 0))
        True

        Different seeds produce different matrices:

        >>> cov1 = rand_cov(3, seed=1)
        >>> cov2 = rand_cov(3, seed=2)
        >>> bool(not np.allclose(cov1, cov2))
        True

        Without seed, consecutive calls may differ (random state):

        >>> # These may or may not be equal depending on random state
        >>> cov_a = rand_cov(2, seed=None)
        >>> cov_b = rand_cov(2, seed=None)
        >>> cov_a.shape == cov_b.shape == (2, 2)
        True

    Note:
        The generated matrix is guaranteed to be positive semi-definite because
        it is constructed as A^T @ A. In practice, it will typically be positive
        definite (all eigenvalues strictly positive) unless n is very large.

    """
    rng = np.random.default_rng(seed)
    a = rng.standard_normal((n, n))
    return np.transpose(a) @ a