qten.pointgroups.ops

Module reference for qten.pointgroups.ops.

ops

Point-group operations on symbolic bases and tensors.

This module combines point-group transforms with QTen Hilbert spaces and tensors. The helpers assemble \(D(g)\) (including the \(SU(2)\) factor on a spinful space), twirl operators by conjugation, and project columns into abelian phase sectors, ordinary finite-group irreps, or projective spinor irreps.

Repository usage

Use hilbert_repr() and the related projection helpers when an existing PointGroupElement, PointGroupOpr, or FinitePointGroup should act on symbolic Hilbert-space data. Group definitions live in qten.pointgroups.elements and qten.pointgroups.finite.

point_group_operator_symmetrize

point_group_operator_symmetrize(
    group: FinitePointGroup,
    operator: Tensor,
    *,
    fixpoint: Offset | None = None,
    rebase_fixpoint: bool = False,
) -> Tensor

Average an operator over a finite point group by unitary conjugation.

This computes [ A_G = |G|^{-1}\sum_{g\in G} D(g)\,A\,D(g)^\dagger. ] The average commutes with every \(D(h)\). Unlike character projection of state vectors, conjugation averaging does not need ordinary or spinor characters: if a lift is replaced by \(\eta(g)u(g)\) with \(\eta(g)\in\{\pm 1\}\), the same factor appears in \(D(g)\) and \(D(g)^\dagger\) and cancels.

The two matrix dimensions of operator must describe the same HilbertSpace. That space must be closed under every point-group operation; otherwise representation assembly raises ValueError.

Parameters:

Name	Type	Description	Default
group	FinitePointGroup	Finite point group whose elements generate the average.	required
operator	Tensor	Rank-2 tensor whose two dimensions are the same ordered HilbertSpace.	required
fixpoint	Offset | None	Desired invariant point. When set, each group element is wrapped as a PointGroupOpr and recentered with fixpoint_at before D(g) is assembled.	None
rebase_fixpoint	bool	Forwarded to PointGroupOpr.fixpoint_at as rebase.	False

Returns:

Type	Description
Tensor	The conjugation-averaged operator on the same Hilbert space.

Raises:

Type	Description
ValueError	If operator is not a square Hilbert-space tensor, or if some D(g) is not numerically unitary.

Source code in src/qten/pointgroups/ops.py

def point_group_operator_symmetrize(
    group: FinitePointGroup,
    operator: Tensor,
    *,
    fixpoint: Offset | None = None,
    rebase_fixpoint: bool = False,
) -> Tensor:
    r"""
    Average an operator over a finite point group by unitary conjugation.

    This computes
    \[
    A_G = |G|^{-1}\sum_{g\in G} D(g)\,A\,D(g)^\dagger.
    \]
    The average commutes with every \(D(h)\). Unlike character projection of
    state vectors, conjugation averaging does not need ordinary or spinor
    characters: if a lift is replaced by \(\eta(g)u(g)\) with
    \(\eta(g)\in\{\pm 1\}\), the same factor appears in \(D(g)\) and
    \(D(g)^\dagger\) and cancels.

    The two matrix dimensions of `operator` must describe the same
    [`HilbertSpace`][qten.symbolics.hilbert_space.HilbertSpace]. That space
    must be closed under every point-group operation; otherwise representation
    assembly raises `ValueError`.

    Parameters
    ----------
    group : FinitePointGroup
        Finite point group whose elements generate the average.
    operator : Tensor
        Rank-2 tensor whose two dimensions are the same ordered
        [`HilbertSpace`][qten.symbolics.hilbert_space.HilbertSpace].
    fixpoint : Offset | None, optional
        Desired invariant point. When set, each group element is wrapped as
        a `PointGroupOpr` and recentered with `fixpoint_at` before `D(g)`
        is assembled.
    rebase_fixpoint : bool, default False
        Forwarded to
        [`PointGroupOpr.fixpoint_at`][qten.pointgroups.elements.PointGroupOpr.fixpoint_at]
        as `rebase`.

    Returns
    -------
    Tensor
        The conjugation-averaged operator on the same Hilbert space.

    Raises
    ------
    ValueError
        If `operator` is not a square Hilbert-space tensor, or if some
        `D(g)` is not numerically unitary.
    """
    if operator.rank() != 2:
        raise ValueError("operator must be a rank-2 Hilbert-space tensor.")
    row_dim, col_dim = operator.dims
    if not isinstance(row_dim, HilbertSpace) or not isinstance(col_dim, HilbertSpace):
        raise ValueError("operator dimensions must both be HilbertSpace objects.")
    if row_dim != col_dim:
        raise ValueError(
            "operator dimensions must use the same ordered Hilbert-space basis, "
            "including identical U(1) gauges."
        )

    elements = group.elements()
    if not elements:
        raise ValueError("Finite point group contains no elements.")

    averaged: Tensor | None = None
    for element in elements:
        element_opr = PointGroupOpr(element)
        if fixpoint is not None:
            element_opr = element_opr.fixpoint_at(fixpoint, rebase=rebase_fixpoint)
        representation = _hilbert_opr_repr(
            element_opr, row_dim, device=operator.device
        ).to_device(operator.device)
        identity = torch.eye(
            representation.data.shape[0],
            dtype=representation.data.dtype,
            device=representation.data.device,
        )
        tolerance = _relative_tolerance(
            representation.data.dtype, representation.data.shape[0]
        )
        gram = representation.data.conj().T @ representation.data
        if not torch.allclose(gram, identity, rtol=tolerance, atol=tolerance):
            max_error = float(torch.max(torch.abs(gram - identity)).item())
            raise ValueError(
                "Point-group operator averaging requires a unitary Hilbert-space "
                f"representation; max |D†D-I|={max_error:.3e} for {element}."
            )
        transformed = representation @ operator @ representation.h(-2, -1)
        averaged = transformed if averaged is None else averaged + transformed

    assert averaged is not None
    return averaged / len(elements)

spinful_transform_basis

spinful_transform_basis(
    opr: PointGroupOpr, psi: U1Basis
) -> U1Span

Apply \(D_{\mathrm{orb}}(g)\otimes u(g)\) to a single basis state.

On a product label \(|\mathrm{orb},s\rangle\), [ D(g)|\mathrm{orb},s\rangle =\sum_{s'}u(g)_{s's}\,|g\cdot\mathrm{orb},\,s'\rangle. ] Spatial irreps that opr allows are transformed as usual. The Spin irrep is expanded with the \(SU(2)\) factor. Returns a U1Span because spin mixing generally produces a superposition.

Source code in src/qten/pointgroups/ops.py

def spinful_transform_basis(opr: PointGroupOpr, psi: U1Basis) -> U1Span:
    r"""
    Apply \(D_{\mathrm{orb}}(g)\otimes u(g)\) to a single basis state.

    On a product label \(|\mathrm{orb},s\rangle\),
    \[
    D(g)|\mathrm{orb},s\rangle
    =\sum_{s'}u(g)_{s's}\,|g\cdot\mathrm{orb},\,s'\rangle.
    \]
    Spatial irreps that `opr` allows are transformed as usual. The
    [`Spin`][qten.phys.spin.Spin] irrep is expanded with the \(SU(2)\) factor.
    Returns a [`U1Span`][qten.symbolics.hilbert_space.U1Span] because spin
    mixing generally produces a superposition.
    """
    nons, spin, psi_coef = _split_spin_basis(psi)
    orbital_terms = _transform_nons_spin_irreps(opr, nons)
    # Standalone transform keeps the raw geometric image. The Hilbert-space
    # representation can fold unit-cell labels for local/Γ-point matching;
    # nonzero-momentum Bloch translation phases are not represented.

    terms: list[U1Basis] = []
    for spatial_coef, new_nons in orbital_terms:
        for amp, spin_out in expand_spin(opr, spin):
            coef = sy.simplify(psi_coef * spatial_coef * amp)
            if coef == 0:
                continue
            terms.append(U1Basis(coef, new_nons + (spin_out,)))

    if not terms:
        raise RuntimeError(f"spinful image of {psi} vanished")

    merged: dict[U1Basis, U1Basis] = {}
    for term in terms:
        ray = term.rays()
        if ray in merged:
            prev = merged[ray]
            merged[ray] = U1Basis(sy.simplify(prev.coef + term.coef), term.base)
        else:
            merged[ray] = term
    return U1Span(tuple(merged.values()))

spinful_hilbert_opr_repr

spinful_hilbert_opr_repr(
    opr: PointGroupOpr,
    space: HilbertSpace,
    *,
    device: Optional[Device] = None,
) -> Tensor

Matrix of \(D(g)=D_{\mathrm{orb}}(g)\otimes u(g)\) on a spinful Hilbert space.

Columns of \(u(g)\) are ordered \((\uparrow,\downarrow)\). Fast path: compute the \(SU(2)\) factor once, cache the orbital image of each distinct non-spin irrep tuple, and scatter numerical amplitudes into a dense matrix. This avoids per-basis SymPy Gram assembly used by the earlier prototype and is suitable for full finite-group symmetrization.

When every Offset label in space is lattice-backed and already intra-cell fractional, orbital images are folded with Offset.fractional() before lookup so primitive-lattice translations match the stored unit-cell basis. This is a local/Γ-point representation; nonzero-momentum Bloch phases require an explicit momentum-dependent representation and are not inserted here.

Parameters:

Name	Type	Description	Default
opr	PointGroupOpr	Point operation, including any affine center already set on it.	required
space	HilbertSpace	Spinful Hilbert space. Every basis state must carry exactly one Spin label, and the space must be closed under opr.	required
device	Optional[Device]	Device for the returned tensor.	None

Returns:

Type	Description
Tensor	Square tensor of \(D_{\mathrm{orb}}(g)\otimes u(g)\) on space.

Source code in src/qten/pointgroups/ops.py

def spinful_hilbert_opr_repr(
    opr: PointGroupOpr,
    space: HilbertSpace,
    *,
    device: Optional[Device] = None,
) -> Tensor:
    r"""
    Matrix of \(D(g)=D_{\mathrm{orb}}(g)\otimes u(g)\) on a spinful Hilbert space.

    Columns of \(u(g)\) are ordered \((\uparrow,\downarrow)\). Fast path:
    compute the \(SU(2)\) factor once, cache the orbital image of each
    distinct non-spin irrep tuple, and scatter numerical amplitudes into a dense
    matrix. This avoids per-basis SymPy Gram assembly used by the earlier
    prototype and is suitable for full finite-group symmetrization.

    When every [`Offset`][qten.geometries.spatials.Offset] label in `space` is
    lattice-backed and already intra-cell fractional, orbital images are
    folded with `Offset.fractional()` before lookup so primitive-lattice
    translations match the stored unit-cell basis. This is a local/Γ-point
    representation; nonzero-momentum Bloch phases require an explicit
    momentum-dependent representation and are not inserted here.

    Parameters
    ----------
    opr : PointGroupOpr
        Point operation, including any affine center already set on it.
    space : HilbertSpace
        Spinful Hilbert space. Every basis state must carry exactly one
        [`Spin`][qten.phys.spin.Spin] label, and the space must be closed
        under `opr`.
    device : Optional[Device], optional
        Device for the returned tensor.

    Returns
    -------
    Tensor
        Square tensor of \(D_{\mathrm{orb}}(g)\otimes u(g)\) on `space`.
    """
    _validate_spinful_space(space)
    precision = get_precision_config()
    torch_device = device.torch_device() if device is not None else None
    n = space.dim
    data = torch.zeros(
        (n, n),
        dtype=precision.torch_complex,
        device=torch_device,
    )

    # u[out, in] in the ordered basis (Spin.up, Spin.down)
    u = torch.tensor(
        su2_numeric(opr),
        dtype=precision.torch_complex,
        device=torch_device,
    )
    spin_list = (Spin.up, Spin.down)
    spin_index = {Spin.up: 0, Spin.down: 1}

    elements = list(space.elements())
    fold_offsets = _space_uses_fractional_offsets(space)
    # Lookup: (non-spin irrep tuple, Spin) -> (basis index, basis coefficient)
    index_of: dict[tuple[tuple[Any, ...], Spin], tuple[int, complex]] = {}
    parsed: list[tuple[tuple[Any, ...], Spin, complex, int]] = []
    for psi in elements:
        nons, spin, coef = _split_spin_basis(psi)
        j = space.structure[psi]
        numeric_coef = complex(sy.N(coef))
        index_of[(nons, spin)] = (j, numeric_coef)
        parsed.append((nons, spin, numeric_coef, j))

    # Cache orbital images of unique non-spin labels
    nons_image: dict[tuple[Any, ...], list[tuple[complex, tuple[Any, ...]]]] = {}
    for nons, _, _, _ in parsed:
        if nons in nons_image:
            continue
        spatial_terms = _transform_nons_spin_irreps(
            opr, nons, space=space, fold_offsets=fold_offsets
        )
        nons_image[nons] = [
            (complex(sy.N(spatial_coef)), new_nons)
            for spatial_coef, new_nons in spatial_terms
        ]

    for nons, spin_in, psi_coef, j in parsed:
        s_in = spin_index[spin_in]
        for spatial_coef, new_nons in nons_image[nons]:
            for s_out, spin_out in enumerate(spin_list):
                amp = u[s_out, s_in]
                if amp.abs().item() == 0.0:
                    continue
                key = (new_nons, spin_out)
                if key not in index_of:
                    raise ValueError(
                        f"Spinful image ({new_nons}, {spin_out}) is not in space "
                        f"{space}. Space is not closed under {opr}."
                    )
                i, target_coef = index_of[key]
                data[i, j] += (target_coef.conjugate() * spatial_coef * psi_coef) * amp

    return Tensor(data=data, dims=(space, space))

get_direct_transform

get_direct_transform(
    opr: PointGroupOpr,
    space: HilbertSpace,
    *,
    device: Optional[Device] = None,
) -> Tensor

Build the external basis-mapping tensor from a Hilbert space to its transformed image.

Unlike hilbert_opr_repr(), this helper does not require opr to preserve the ray structure of space. Instead it explicitly constructs the transformed output HilbertSpace and returns a one-hot mapping matrix with dims (space, out_space).

Spinful Hilbert spaces are rejected: a generic \(SU(2)\) factor maps one basis state to a superposition, which this one-to-one mapping cannot express. Use hilbert_repr instead.

When a basis state contains a PointGroupBasis irrep, that irrep is transformed directly in the Euclidean polynomial basis. In particular, no eigen-phase is factored out. For example, a basis function x rotated by C4 is mapped to y in the output space rather than left as x with a phase in the tensor data.

Parameters:

Name	Type	Description	Default
opr	PointGroupOpr	Point-group operator used to transform basis labels.	required
space	HilbertSpace	Input Hilbert space whose ordered basis defines the source axis.	required
device	Optional[Device]	Device on which to allocate the returned mapping tensor.	None

Returns:

Type	Description
Tensor	Rank-2 tensor with dimensions (space, out_space) and only 1 numerical entries at the mapped basis positions.

Raises:

Type	Description
NotImplementedError	If any basis state in space carries a spin-1/2 label.

Source code in src/qten/pointgroups/ops.py

def get_direct_transform(
    opr: PointGroupOpr,
    space: HilbertSpace,
    *,
    device: Optional[Device] = None,
) -> Tensor:
    r"""
    Build the external basis-mapping tensor from a Hilbert space to its transformed image.

    Unlike [`hilbert_opr_repr()`][qten.symbolics.ops.hilbert_opr_repr], this helper does not require `opr` to preserve the ray structure of
    `space`. Instead it explicitly constructs the transformed output
    [`HilbertSpace`][qten.symbolics.hilbert_space.HilbertSpace] and returns a one-hot mapping matrix with dims `(space, out_space)`.

    Spinful Hilbert spaces are rejected: a generic \(SU(2)\) factor maps one
    basis state to a superposition, which this one-to-one mapping cannot
    express. Use [`hilbert_repr`][qten.pointgroups.ops.hilbert_repr] instead.

    When a basis state contains a [`PointGroupBasis`][qten.pointgroups.basis.PointGroupBasis] irrep, that irrep is transformed directly in the Euclidean polynomial basis.
    In particular, no eigen-phase is factored out. For example, a basis
    function `x` rotated by `C4` is mapped to `y` in the output space rather
    than left as `x` with a phase in the tensor data.

    Parameters
    ----------
    opr : PointGroupOpr
        Point-group operator used to transform basis labels.
    space : HilbertSpace
        Input Hilbert space whose ordered basis defines the source axis.
    device : Optional[Device], optional
        Device on which to allocate the returned mapping tensor.

    Returns
    -------
    Tensor
        Rank-2 tensor with dimensions `(space, out_space)` and only `1`
        numerical entries at the mapped basis positions.

    Raises
    ------
    NotImplementedError
        If any basis state in `space` carries a spin-1/2 label.
    """
    transformed = {psi: _ext_transform_basis(opr, psi) for psi in space.elements()}
    out_space = space.map(lambda psi: transformed[psi])
    return mapping_matrix(space, out_space, transformed, device=device)

hilbert_repr

hilbert_repr(
    opr: PointGroupOpr,
    space: HilbertSpace,
    *,
    device: Optional[Device] = None,
) -> Tensor

Assemble the Hilbert-space representation \(D(g)\) of a point operation.

Spinless spaces use \(D(g)=D_{\mathrm{orb}}(g)\). Spinful spaces use [ D(g)\,|\mathrm{orb},s\rangle =\sum_{s'}u(g)_{s's}\,|g\cdot\mathrm{orb},\,s'\rangle, ] i.e. \(D(g)=D_{\mathrm{orb}}(g)\otimes u(g)\). Each basis irrep is transformed on its own: lattice Offset labels may be folded back into the unit cell, PointGroupBasis polynomials are canonicalized onto labels already present in space, and Spin is expanded by the \(SU(2)\) lift. This function has no fixpoint=; recenter opr with fixpoint_at first.

Parameters:

Name	Type	Description	Default
opr	PointGroupOpr	Point operation, including any affine center already set on it.	required
space	HilbertSpace	Ordered basis. Must be closed under opr.	required
device	Optional[Device]	Device for the returned tensor.	None

Returns:

Type	Description
Tensor	Square tensor of \(D(g)\) on space. On a spinful space this is \(D_{\mathrm{orb}}(g)\otimes u(g)\).

Source code in src/qten/pointgroups/ops.py

def hilbert_repr(
    opr: PointGroupOpr, space: HilbertSpace, *, device: Optional[Device] = None
) -> Tensor:
    r"""
    Assemble the Hilbert-space representation \(D(g)\) of a point operation.

    Spinless spaces use \(D(g)=D_{\mathrm{orb}}(g)\). Spinful spaces use
    \[
    D(g)\,|\mathrm{orb},s\rangle
    =\sum_{s'}u(g)_{s's}\,|g\cdot\mathrm{orb},\,s'\rangle,
    \]
    i.e. \(D(g)=D_{\mathrm{orb}}(g)\otimes u(g)\). Each basis irrep is
    transformed on its own: lattice
    [`Offset`][qten.geometries.spatials.Offset] labels may be folded back into
    the unit cell, [`PointGroupBasis`][qten.pointgroups.basis.PointGroupBasis]
    polynomials are canonicalized onto labels already present in `space`, and
    [`Spin`][qten.phys.spin.Spin] is expanded by the \(SU(2)\) lift. This
    function has no `fixpoint=`; recenter `opr` with
    [`fixpoint_at`][qten.pointgroups.elements.PointGroupOpr.fixpoint_at]
    first.

    Parameters
    ----------
    opr : PointGroupOpr
        Point operation, including any affine center already set on it.
    space : HilbertSpace
        Ordered basis. Must be closed under `opr`.
    device : Optional[Device], optional
        Device for the returned tensor.

    Returns
    -------
    Tensor
        Square tensor of \(D(g)\) on `space`. On a spinful space this is
        \(D_{\mathrm{orb}}(g)\otimes u(g)\).
    """
    _validate_hilbert_basis(space)
    if contains_spin(space):
        return spinful_hilbert_opr_repr(opr, space, device=device)
    if not _contains_point_group_basis(space):
        return hilbert_opr_repr(opr, space, device=device)

    ray_to_basis = {psi.rays(): psi for psi in space.elements()}
    precision = get_precision_config()
    torch_device = device.torch_device() if device is not None else None
    data = torch.zeros(
        (space.dim, space.dim),
        dtype=precision.torch_complex,
        device=torch_device,
    )

    fold_offsets = _space_uses_fractional_offsets(space)
    for source in space.elements():
        images = _transform_nons_spin_irreps(
            opr, source.base, space=space, fold_offsets=fold_offsets
        )
        j = space.structure[source]
        for spatial_coef, new_base in images:
            constructed = U1Basis(sy.simplify(source.coef * spatial_coef), new_base)
            target = ray_to_basis.get(constructed.rays())
            if target is None:
                raise ValueError("opr does not preserve the ray structure of space.")
            i = space.structure[target]
            matrix_element = sy.simplify(sy.conjugate(target.coef) * constructed.coef)
            data[i, j] += complex(sy.N(matrix_element))

    return Tensor(data=data, dims=(space, space))

joint_point_group_basis

joint_point_group_basis(
    oprs: Sequence[PointGroupElement | PointGroupOpr],
    order: int,
) -> FrozenDict[
    tuple[sy.Expr, ...], tuple[PointGroupBasis, ...]
]

Compute common Euclidean eigenfunctions for a commuting family of abelian operators.

The returned table is keyed by one phase per input operator. Each value is the tuple of normalized PointGroupBasis functions spanning the simultaneous eigenspace for that joint phase sector.

Parameters:

Name	Type	Description	Default
oprs	Sequence[PointGroupElement | PointGroupOpr]	Non-empty sequence of operators. Affine PointGroupOpr inputs contribute only their linear part.	required
order	int	Homogeneous polynomial degree used for all Euclidean representations.	required

Returns:

Type	Description
FrozenDict[tuple[Expr, ...], tuple[PointGroupBasis, ...]]	Mapping from joint phase tuple to the simultaneous eigen-basis functions for that sector.

Raises:

Type	Description
ValueError	If oprs is empty, if the operators do not share the same ordered axes, or if their Euclidean representations at order do not commute.

Source code in src/qten/pointgroups/ops.py

def joint_point_group_basis(
    oprs: Sequence[PointGroupElement | PointGroupOpr], order: int
) -> FrozenDict[tuple[sy.Expr, ...], tuple[PointGroupBasis, ...]]:
    """
    Compute common Euclidean eigenfunctions for a commuting family of abelian operators.

    The returned table is keyed by one phase per input operator. Each value is
    the tuple of normalized [`PointGroupBasis`][qten.pointgroups.basis.PointGroupBasis]
    functions spanning the simultaneous eigenspace for that joint phase sector.

    Parameters
    ----------
    oprs : Sequence[PointGroupElement | PointGroupOpr]
        Non-empty sequence of operators. Affine
        [`PointGroupOpr`][qten.pointgroups.elements.PointGroupOpr] inputs contribute
        only their linear part.
    order : int
        Homogeneous polynomial degree used for all Euclidean representations.

    Returns
    -------
    FrozenDict[tuple[sy.Expr, ...], tuple[PointGroupBasis, ...]]
        Mapping from joint phase tuple to the simultaneous eigen-basis
        functions for that sector.

    Raises
    ------
    ValueError
        If `oprs` is empty, if the operators do not share the same ordered
        axes, or if their Euclidean representations at `order` do not commute.
    """
    if not oprs:
        raise ValueError("oprs must be non-empty.")

    groups = tuple(opr.g if isinstance(opr, PointGroupOpr) else opr for opr in oprs)
    axes = groups[0].axes
    if any(g.axes != axes for g in groups[1:]):
        raise ValueError("All operators must share the same ordered axes.")

    transforms = tuple(g.euclidean_repr(order) for g in groups)
    zero = sy.zeros(transforms[0].rows, transforms[0].cols)
    for i, left in enumerate(transforms):
        for right in transforms[i + 1 :]:
            if not sy.simplify(left @ right - right @ left).equals(zero):
                raise ValueError(
                    "All operators must commute in the Euclidean representation "
                    f"of order {order}."
                )

    euclidean_basis = groups[0].euclidean_basis(order)
    ident = sy.ImmutableDenseMatrix.eye(transforms[0].rows)
    all_sector_projectors: list[list[tuple[sy.Expr, sy.ImmutableDenseMatrix]]] = []
    for g, transform in zip(groups, transforms):
        powers = [ident]
        for _ in range(1, g.group_order()):
            powers.append(sy.ImmutableDenseMatrix(sy.simplify(powers[-1] @ transform)))

        sector_projectors: list[tuple[sy.Expr, sy.ImmutableDenseMatrix]] = []
        for phase in g.basis(order):
            projector = sy.zeros(transform.rows, transform.cols)
            for k, power in enumerate(powers):
                projector += sy.simplify((phase ** (-k)) * power)
            sector_projectors.append(
                (
                    sy.simplify(phase),
                    sy.ImmutableDenseMatrix(sy.simplify(projector / g.group_order())),
                )
            )
        all_sector_projectors.append(sector_projectors)

    tbl: dict[tuple[sy.Expr, ...], tuple[PointGroupBasis, ...]] = {}
    for sector_product in product(*all_sector_projectors):
        phases = tuple(phase for phase, _ in sector_product)
        projector = ident
        for _, sector_projector in sector_product:
            projector = sy.ImmutableDenseMatrix(
                sy.simplify(sector_projector @ projector)
            )

        basis_vectors = projector.columnspace()
        if not basis_vectors:
            continue

        labels: list[PointGroupBasis] = []
        seen_reps = set()
        for vec in basis_vectors:
            rep = sy.ImmutableDenseMatrix(vec)
            if all(entry == 0 for entry in rep):
                continue
            basis = PointGroupBasis.from_rep(
                rep=rep,
                euclidean_basis=euclidean_basis,
                axes=axes,
                order=order,
                irrep=phases,
            )
            rep_key = tuple(basis.rep)
            if rep_key in seen_reps:
                continue
            seen_reps.add(rep_key)
            labels.append(basis)

        if labels:
            tbl[phases] = tuple(labels)

    return FrozenDict(tbl)

point_group_column_symmetrize

point_group_column_symmetrize(
    opr: PointGroupOpr | FinitePointGroup,
    w: Tensor,
    full_sector: bool = False,
    *,
    fixpoint: Offset | None = None,
    rebase_fixpoint: bool = False,
) -> Tensor

Symmetrize the columns of w by projecting each one onto symmetry sectors.

For a finite-order abelian operator opr of spatial order \(n\), each exact spinless sector is labeled by a root of unity \(\zeta^n=1\). For spin-1/2 the safe common period is \(N=2n\), because a \(2\pi\) proper spin rotation is \(-I\), and the sector is a SpinfulPhaseSector. The projector on \(G=D(g)\) is [ P_\zeta=\frac{1}{N}\sum_{k=0}^{N-1}\zeta^{-k}G^k,\qquad\zeta^N=1, ] with \(N=n\) spinless and \(N=2n\) spinful. The period \(2n\) need not be minimal when the proper spin factor is already the identity.

If opr is a FinitePointGroup, this routine uses [ P^\mu=\frac{d_\mu}{|G|}\sum_{g\in G}\chi^\mu(g)^*D(g). ] Spinless spaces use ordinary (linear) \(\chi\). Spaces that already contain Spin use the group's \(SU(2)\) section and element-wise projective \(\chi\), unless the group was defined with spin="trivial". A cyclic PointGroupOpr is the abelian special case: one-dimensional characters \(\zeta^k\).

The projector is applied to each input column separately. When full_sector is True, every nonzero projected sector component is returned. When full_sector is False, only the dominant nonzero sector component of each input column is kept, so the output column count does not exceed the input count. Returned columns carry a sector label: FiniteIrrepSector for ordinary finite-group irreps, SpinorIrrepSector for projective spinor irreps, or SpinfulPhaseSector / PointGroupBasis for abelian phase sectors.

The output column count can differ from the input one only when full_sector=True, because symmetry projection may split one approximate column into multiple exact sectors.

Parameters:

Name	Type	Description	Default
opr	PointGroupOpr | FinitePointGroup	Symmetry descriptor. PointGroupOpr uses the abelian phase-sector path; FinitePointGroup uses finite-group irrep projectors.	required
w	Tensor	Rank-2 tensor whose first dimension is a HilbertSpace and whose columns are vectors to project.	required
full_sector	bool	If True, return every nonzero sector component of each input column. If False, keep only the largest nonzero sector component per input column.	False
fixpoint	Offset | None	Desired invariant point. When set, each group element is wrapped as a PointGroupOpr and recentered with fixpoint_at before D(g) is assembled.	None
rebase_fixpoint	bool	Forwarded to PointGroupOpr.fixpoint_at as rebase.	False

Returns:

Type	Description
Tensor	Rank-2 tensor with the same row Hilbert space and a column HilbertSpace labelled by symmetry-sector basis data.

Raises:

Type	Description
ValueError	If w is not rank 2, if w.dims[0] is not a HilbertSpace, or if w.dims[1] is neither an IndexSpace nor a HilbertSpace.

Source code in src/qten/pointgroups/ops.py

def point_group_column_symmetrize(
    opr: PointGroupOpr | FinitePointGroup,
    w: Tensor,
    full_sector: bool = False,
    *,
    fixpoint: Offset | None = None,
    rebase_fixpoint: bool = False,
) -> Tensor:
    r"""
    Symmetrize the columns of `w` by projecting each one onto symmetry sectors.

    For a finite-order abelian operator `opr` of spatial order \(n\), each
    exact spinless sector is labeled by a root of unity \(\zeta^n=1\). For
    spin-1/2 the safe common period is \(N=2n\), because a \(2\pi\) proper
    spin rotation is \(-I\), and the sector is a
    [`SpinfulPhaseSector`][qten.pointgroups.sectors.SpinfulPhaseSector].
    The projector on \(G=D(g)\) is
    \[
    P_\zeta=\frac{1}{N}\sum_{k=0}^{N-1}\zeta^{-k}G^k,\qquad\zeta^N=1,
    \]
    with \(N=n\) spinless and \(N=2n\) spinful. The period \(2n\) need not
    be minimal when the proper spin factor is already the identity.

    If `opr` is a [`FinitePointGroup`][qten.pointgroups.finite.FinitePointGroup],
    this routine uses
    \[
    P^\mu=\frac{d_\mu}{|G|}\sum_{g\in G}\chi^\mu(g)^*D(g).
    \]
    Spinless spaces use ordinary (linear) \(\chi\). Spaces that already
    contain [`Spin`][qten.phys.spin.Spin] use the group's \(SU(2)\) section and
    element-wise projective \(\chi\), unless the group was defined with
    `spin="trivial"`. A cyclic
    [`PointGroupOpr`][qten.pointgroups.elements.PointGroupOpr] is the
    abelian special case: one-dimensional characters \(\zeta^k\).

    The projector is applied to each input column separately. When
    `full_sector` is `True`, every
    nonzero projected sector component is returned. When `full_sector` is
    `False`, only the dominant nonzero sector component of each input column is
    kept, so the output column count does not exceed the input count.
    Returned columns carry a sector label:
    [`FiniteIrrepSector`][qten.pointgroups.sectors.FiniteIrrepSector] for ordinary
    finite-group irreps,
    [`SpinorIrrepSector`][qten.pointgroups.sectors.SpinorIrrepSector] for
    projective spinor irreps, or
    [`SpinfulPhaseSector`][qten.pointgroups.sectors.SpinfulPhaseSector] /
    [`PointGroupBasis`][qten.pointgroups.basis.PointGroupBasis] for abelian
    phase sectors.

    The output column count can differ from the input one only when
    `full_sector=True`, because symmetry projection may split one approximate
    column into multiple exact sectors.

    Parameters
    ----------
    opr : PointGroupOpr | FinitePointGroup
        Symmetry descriptor. `PointGroupOpr` uses the abelian phase-sector path;
        `FinitePointGroup` uses finite-group irrep projectors.
    w : Tensor
        Rank-2 tensor whose first dimension is a
        [`HilbertSpace`][qten.symbolics.hilbert_space.HilbertSpace] and whose
        columns are vectors to project.
    full_sector : bool, default False
        If `True`, return every nonzero sector component of each input column.
        If `False`, keep only the largest nonzero sector component per input
        column.
    fixpoint : Offset | None, optional
        Desired invariant point. When set, each group element is wrapped as
        a `PointGroupOpr` and recentered with `fixpoint_at` before `D(g)`
        is assembled.
    rebase_fixpoint : bool, default False
        Forwarded to
        [`PointGroupOpr.fixpoint_at`][qten.pointgroups.elements.PointGroupOpr.fixpoint_at]
        as `rebase`.

    Returns
    -------
    Tensor
        Rank-2 tensor with the same row Hilbert space and a column
        [`HilbertSpace`][qten.symbolics.hilbert_space.HilbertSpace] labelled by
        symmetry-sector basis data.

    Raises
    ------
    ValueError
        If `w` is not rank 2, if `w.dims[0]` is not a `HilbertSpace`, or if
        `w.dims[1]` is neither an `IndexSpace` nor a `HilbertSpace`.
    """
    if isinstance(opr, FinitePointGroup):
        return _finite_point_group_column_symmetrize(
            opr,
            w,
            full_sector=full_sector,
            fixpoint=fixpoint,
            rebase_fixpoint=rebase_fixpoint,
        )

    if fixpoint is not None:
        opr = opr.fixpoint_at(fixpoint, rebase=rebase_fixpoint)

    row_dim, seeds = _column_symmetrize_context(w)

    g_full = _hilbert_opr_repr(opr, row_dim, device=w.device).to_device(w.device)
    spatial_order = opr.g.group_order()
    spinful = contains_spin(row_dim)
    order = 2 * spatial_order if spinful else spatial_order
    ident = eye((row_dim, row_dim)).astype(g_full.data.dtype).to_device(g_full.device)
    single_col = IndexSpace.linear(1)
    tol = _relative_tolerance(g_full.data.dtype, row_dim.dim)

    g_powers: list[Tensor] = [ident]
    for _ in range(1, order):
        g_powers.append(g_powers[-1] @ g_full)

    closure = g_powers[-1] @ g_full
    if not torch.allclose(closure.data, ident.data, rtol=0.0, atol=tol):
        raise ValueError(
            f"Operator representation does not close at expected order {order}"
        )

    sector_projectors: list[tuple[Any, Tensor]] = []
    for m in range(order):
        phase_exact = sy.simplify(sy.exp(2 * sy.pi * sy.I * m / order))
        sector_label: Any
        if spinful:
            sector_label = SpinfulPhaseSector(
                phase=phase_exact,
                spatial_order=spatial_order,
            )
        else:
            sector_label = _phase_basis(opr, phase_exact)
        phase_scalar = complex(sy.N(phase_exact))

        projector = 0 * ident
        for k, g_power in enumerate(g_powers):
            projector = projector + (phase_scalar ** (-k)) * g_power
        sector_projectors.append((sector_label, projector / order))

    projected_cols: list[Tensor] = []
    raw_labels: list[U1Basis] = []
    for j, seed in enumerate(seeds):
        col = (
            w[:, j : j + 1].clone().replace_dim(1, single_col).astype(g_full.data.dtype)
        )
        input_norm = abs(col.norm().item())
        projection_cutoff = tol * input_norm
        candidates: list[tuple[float, Tensor, U1Basis]] = []
        for sector_label, projector in sector_projectors:
            projected = projector @ col

            projected_norm = projected.norm()
            norm_value = abs(projected_norm.item())
            if norm_value <= projection_cutoff:
                continue

            candidates.append(
                (
                    norm_value,
                    projected / norm_value,
                    _attach_sector_label(seed, sector_label),
                )
            )

        if full_sector:
            for _, projected, label in candidates:
                projected_cols.append(projected)
                raw_labels.append(label)
        elif candidates:
            _, projected, label = max(candidates, key=lambda item: item[0])
            projected_cols.append(projected)
            raw_labels.append(label)

    if not projected_cols:
        return Tensor(
            data=w.data.new_empty((row_dim.dim, 0), dtype=g_full.data.dtype),
            dims=(row_dim, IndexSpace.linear(0)),
        )

    out_dim = HilbertSpace.new(_labels_with_degeneracy(raw_labels))
    return cat(projected_cols, dim=-1).replace_dim(-1, out_dim)

joint_point_group_column_symmetrize

joint_point_group_column_symmetrize(
    oprs: Sequence[PointGroupOpr],
    w: Tensor,
    full_sector: bool = False,
    *,
    group: FinitePointGroup | None = None,
) -> Tensor

Symmetrize columns of w into simultaneous sectors of abelian operators.

The operators in oprs are expected to commute on w.dims[0]. For each operator this builds the same \(P_\zeta\) as point_group_column_symmetrize. A joint sector is the product projector \(P_{\zeta_1}\cdots P_{\zeta_m}\) over the Cartesian product of those roots of unity.

When full_sector is True, every nonzero joint-sector component is returned. When False, only the dominant nonzero joint-sector component of each input column is kept. Spinless columns carry a representative common PointGroupBasis for the corresponding joint phase sector. Spinful columns carry one JointSpinfulPhaseSector.

This helper has no fixpoint= argument. Center each operator with fixpoint_at before calling it.

Parameters:

Name	Type	Description	Default
oprs	Sequence[PointGroupOpr]	Non-empty sequence of finite-order abelian operators. They are expected to commute on the row Hilbert space of w.	required
w	Tensor	Rank-2 tensor whose first dimension is a HilbertSpace and whose columns are vectors to project.	required
full_sector	bool	If True, return every nonzero joint-sector component of each input column. If False, keep only the largest nonzero joint-sector component per input column.	False
group	FinitePointGroup | None	Parent finite point group. Required when w is spinful: every operator must be an element of this group, so the joint family stays inside one double cover.	None

Returns:

Type	Description
Tensor	Rank-2 tensor with the same row Hilbert space and a column HilbertSpace labelled by representative joint-sector basis data.

Raises:

Type	Description
ValueError	If oprs is empty, if w is not rank 2, if w.dims[0] is not a HilbertSpace, or if w.dims[1] is neither an IndexSpace nor a HilbertSpace.

Source code in src/qten/pointgroups/ops.py

def joint_point_group_column_symmetrize(
    oprs: Sequence[PointGroupOpr],
    w: Tensor,
    full_sector: bool = False,
    *,
    group: FinitePointGroup | None = None,
) -> Tensor:
    r"""
    Symmetrize columns of `w` into simultaneous sectors of abelian operators.

    The operators in `oprs` are expected to commute on `w.dims[0]`. For each
    operator this builds the same \(P_\zeta\) as
    [`point_group_column_symmetrize`][qten.pointgroups.ops.point_group_column_symmetrize].
    A joint sector is the product projector
    \(P_{\zeta_1}\cdots P_{\zeta_m}\) over the Cartesian product of those
    roots of unity.

    When `full_sector` is `True`, every nonzero joint-sector component is
    returned. When `False`, only the dominant nonzero joint-sector component of
    each input column is kept. Spinless columns carry a representative common
    [`PointGroupBasis`][qten.pointgroups.basis.PointGroupBasis] for the
    corresponding joint phase sector. Spinful columns carry one
    [`JointSpinfulPhaseSector`][qten.pointgroups.sectors.JointSpinfulPhaseSector].

    This helper has no `fixpoint=` argument. Center each operator with
    [`fixpoint_at`][qten.pointgroups.elements.PointGroupOpr.fixpoint_at]
    before calling it.

    Parameters
    ----------
    oprs : Sequence[PointGroupOpr]
        Non-empty sequence of finite-order abelian operators. They are expected
        to commute on the row Hilbert space of `w`.
    w : Tensor
        Rank-2 tensor whose first dimension is a
        [`HilbertSpace`][qten.symbolics.hilbert_space.HilbertSpace] and whose
        columns are vectors to project.
    full_sector : bool, default False
        If `True`, return every nonzero joint-sector component of each input
        column. If `False`, keep only the largest nonzero joint-sector
        component per input column.
    group : FinitePointGroup | None, optional
        Parent finite point group. Required when `w` is spinful: every operator
        must be an element of this group, so the joint family stays inside one
        double cover.

    Returns
    -------
    Tensor
        Rank-2 tensor with the same row Hilbert space and a column
        [`HilbertSpace`][qten.symbolics.hilbert_space.HilbertSpace] labelled by
        representative joint-sector basis data.

    Raises
    ------
    ValueError
        If `oprs` is empty, if `w` is not rank 2, if `w.dims[0]` is not a
        `HilbertSpace`, or if `w.dims[1]` is neither an `IndexSpace` nor a
        `HilbertSpace`.
    """
    if not oprs:
        raise ValueError("oprs must be non-empty.")
    if len(oprs) == 1:
        return point_group_column_symmetrize(oprs[0], w, full_sector=full_sector)
    row_dim, seeds = _column_symmetrize_context(w)
    spinful = contains_spin(row_dim)
    if group is not None:
        _require_joint_oprs_in_group(oprs, group)
    elif spinful:
        raise ValueError(
            "Joint spinful projection requires group= the FinitePointGroup "
            "that contains every operator."
        )

    single_col = IndexSpace.linear(1)
    joint_sector_bases = None if spinful else _joint_phase_basis(oprs)
    all_sector_projectors: list[list[tuple[sy.Expr, Tensor]]] = []
    representations: list[Tensor] = []
    dtype = w.data.dtype
    device = w.device
    for opr in oprs:
        g_full = _hilbert_opr_repr(opr, row_dim, device=w.device).to_device(w.device)
        dtype = g_full.data.dtype
        device = g_full.device
        tolerance = _relative_tolerance(dtype, row_dim.dim)
        for previous in representations:
            commutator = g_full.data @ previous.data - previous.data @ g_full.data
            if not torch.allclose(
                commutator,
                torch.zeros_like(commutator),
                rtol=0.0,
                atol=tolerance,
            ):
                max_error = float(torch.max(torch.abs(commutator)).item())
                raise ValueError(
                    "Joint point-group projection requires commuting Hilbert-space "
                    f"representations; max |[D_i,D_j]|={max_error:.3e}."
                )
        representations.append(g_full)
        spatial_order = opr.g.group_order()
        order = 2 * spatial_order if spinful else spatial_order
        ident = eye((row_dim, row_dim)).astype(dtype).to_device(device)

        g_powers: list[Tensor] = [ident]
        for _ in range(1, order):
            g_powers.append(g_powers[-1] @ g_full)
        closure = g_powers[-1] @ g_full
        if not torch.allclose(closure.data, ident.data, rtol=0.0, atol=tolerance):
            raise ValueError(
                f"Operator representation does not close at expected order {order}"
            )

        sector_projectors: list[tuple[sy.Expr, Tensor]] = []
        for m in range(order):
            phase_exact = sy.simplify(sy.exp(2 * sy.pi * sy.I * m / order))
            phase_scalar = complex(sy.N(phase_exact))

            projector = 0 * ident
            for k, g_power in enumerate(g_powers):
                projector = projector + (phase_scalar ** (-k)) * g_power
            sector_projectors.append((phase_exact, projector / order))
        all_sector_projectors.append(sector_projectors)

    projected_cols: list[Tensor] = []
    raw_labels: list[U1Basis] = []
    for j, seed in enumerate(seeds):
        col = w[:, j : j + 1].clone().replace_dim(1, single_col).astype(dtype)
        input_norm = abs(col.norm().item())
        projection_cutoff = _relative_tolerance(dtype, row_dim.dim) * input_norm
        candidates: list[tuple[float, Tensor, U1Basis]] = []
        for sector_product in product(*all_sector_projectors):
            phases = tuple(sy.simplify(phase) for phase, _ in sector_product)
            if spinful:
                label = _attach_sector_label(
                    seed,
                    JointSpinfulPhaseSector(
                        phases=phases,
                        spatial_orders=tuple(opr.g.group_order() for opr in oprs),
                    ),
                )
            else:
                assert joint_sector_bases is not None
                basis = joint_sector_bases.get(phases)
                if basis is None:
                    continue
                label = _attach_basis_label(seed, basis)
            projected = col
            for _, projector in sector_product:
                projected = projector @ projected

            projected_norm = projected.norm()
            norm_value = abs(projected_norm.item())
            if norm_value <= projection_cutoff:
                continue

            candidates.append((norm_value, projected / norm_value, label))

        if full_sector:
            for _, projected, label in candidates:
                projected_cols.append(projected)
                raw_labels.append(label)
        elif candidates:
            _, projected, label = max(candidates, key=lambda item: item[0])
            projected_cols.append(projected)
            raw_labels.append(label)

    if not projected_cols:
        return Tensor(
            data=w.data.new_empty((row_dim.dim, 0), dtype=dtype),
            dims=(row_dim, IndexSpace.linear(0)),
        )

    out_dim = HilbertSpace.new(_labels_with_degeneracy(raw_labels))
    return cat(projected_cols, dim=-1).replace_dim(-1, out_dim)