# Lesson 16 — Loading a useful trial state

## The state we want and the state we can afford

We have an active Hamiltonian for stretched nitrogen: neutral singlet N₂ at 1.85 Å, reduced from the original 28-spatial-orbital cc-pVDZ model to a six-electron, six-spatial-orbital model. Its ground eigenstate is the scientific target, but knowing the operator does not mean knowing that eigenstate’s amplitudes.

A **trial state** is a state we can construct with available classical information and a practical circuit. A Hartree–Fock determinant is cheap but may poorly describe stretched bonds. A selected-CI state can improve the overlap without solving the full original problem. A complete CASCI state in a small active space is a useful reference, but an exact classical solution of every large proposed target is not a scalable prerequisite for quantum computation.

State preparation has two costs: obtaining the classical description and implementing that description coherently. Four hundred nonzero CI coefficients are not four hundred free physical instructions, and a list of amplitudes in memory is not a quantum state already loaded into a processor.

This lesson traces every current registered `state_prep` family, plus the identity helper and the sparse method’s optional binary compression. It also separates chemical trial-state loading from coefficient-index preparation for block encodings. These routines share some software interfaces but do not all produce the same kind of output state. All examples are unexecuted prose; existing source tests are evidence of intended behavior, not tests run here.

## Ansatz, coefficients, and the alignment contract

An **ansatz** usually means a chosen form of approximate state. In this repository, [`data.Ansatz`](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/include/qdk/chemistry/data/ansatz.hpp) has a concrete narrower meaning: it bundles a `Hamiltonian` with a `Wavefunction`. It is not necessarily a parameterized quantum circuit or a variational optimizer.

The [implementation](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/data/ansatz.cpp) validates the pairing. If both objects share the same orbital object, it accepts that identity directly. Otherwise it compares MO/AO dimensions, restricted/unrestricted status, alpha/beta active and inactive indices, orbital coefficient matrices, and available orbital energies. Coefficient-matrix differences are compared against machine epsilon in the norm used there; available energies use `1e-12`. Two bases spanning mathematically equivalent spaces but differing by rotations or phase choices are not automatically reconciled. This is deliberate alignment checking, not automatic orbital transport.

`Ansatz.calculate_energy()` is implemented despite a stale header note suggesting future implementation. It requires available RDMs and contracts them with integrals, adding the core energy. The current implementation uses real RDM variants and has specific restricted/unrestricted block contracts. It cannot turn arbitrary stored CC amplitudes into an expectation value by wishful thinking. The projected-MC eager spin-traced two-RDM normalization caveat from Lesson 9 also still matters: do not use that particular payload as an unchecked energy-contraction reference.

For state preparation, the relevant description is normally

\[
|\psi\rangle=\sum_{j\in\mathcal S}c_j|D_j\rangle,
\qquad \sum_j|c_j|^2=1.
\]

Here \(D_j\) is a determinant in one specified orbital/mode ordering, \(c_j\) its complex amplitude, and \(\mathcal S\) its **support**, the set of nonzero components. **Sparsity** means that support is small relative to all occupation bitstrings. Its usefulness depends not only on the number of determinants but also on the relationships among their bits.

A signed coefficient is not a probability. For real amplitudes, minus signs affect interference. For complex amplitudes, relative phases matter too. Replacing all coefficients with their absolute values generally changes the state, even though direct occupation measurements give the same probabilities.

Only containers exposing actual determinant coefficients can feed these wavefunction loaders. `AmplitudeContainer` stores cluster amplitudes, not a ready determinant expansion; localizers and active-space selectors often return determinant carriers rather than the original correlated coefficients. Obtain a compatible correlated state in the selected, possibly transformed orbitals before treating it as the trial CI state.

## Overlap, truncation, and reoptimization

For normalized states, **fidelity** means squared overlap,

\[
F=|\langle\psi_{\mathrm{target}}|\psi_{\mathrm{trial}}\rangle|^2.
\]

Some source tests name the unsquared overlap magnitude “fidelity”; we will keep the squared definition and identify what a test actually checks. Overlap requires consistent orbital basis, determinant identities, signs, and encoding. The wavefunction container’s overlap routine has its own ordering/pointer restrictions; equal vector lengths are not sufficient evidence of matching determinant bases.

If a normalized expansion is truncated to a chosen support \(T\), its retained probability mass is \(w=\sum_{j\in T}|c_j|^2\). Renormalization divides retained coefficients by \(\sqrt w\). The new state has fidelity \(w\) with the original expansion, provided determinants are orthonormal and alignment is unchanged. This is not automatically fidelity \(w\) with the unknown exact ground state.

The `Wavefunction.truncate` operation selects and renormalizes amplitudes. It does not solve for better amplitudes in the chosen support. A projected multiconfiguration calculation instead forms the Hamiltonian in that selected determinant space and solves its eigenproblem, as taught in Lesson 9.

A hand-worked example separates them. Suppose a selected two-determinant subspace has

\[
H_T=\begin{pmatrix}0&-1\\-1&0\end{pmatrix}.
\]

A truncated-and-renormalized trial vector \((3,1)/\sqrt{10}\) has expectation −0.6 in these illustrative energy units. Reoptimizing the coefficients gives \((1,1)/\sqrt2\) and eigenvalue −1. No determinants were added; their mixing was improved. This is a toy matrix, not a nitrogen result and not an executed PMC calculation.

For a Hermitian model and a normalized state, the variational expectation is bounded below by that model’s ground energy. It need not be an eigenvalue: a state spread over eigenstates has their probability-weighted average energy. Enlarging a nested variational determinant space cannot raise its exactly optimized minimum, but it can increase preprocessing and loading cost, and does not guarantee a better end-to-end resource tradeoff.

If an ideal coherent phase-estimation procedure acts on a trial state \(\sum_k a_k|E_k\rangle\), it identifies the eigenstate branch \(k\) with weight \(|a_k|^2\), subject to its phase resolution. Under that idealized experiment, ground-state branch probability equals ground-state overlap weight. The repository’s iterative bit-estimation circuits can freshly prepare a state for each bit/shot. Their full bit-string success probability is not generically “the fidelity”: fresh eigenstate mixtures, finite precision, feedback, and repeated sampling must be accounted for. We will derive that distinction in the phase-estimation lesson.

## The offered loaders and their real inputs

The [registry](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/algorithms/registry.py) registers `dense_pure_state`, `sparse_isometry`, `qrom`, and `alias_sampling`; the default is `sparse_isometry`. The optional Qiskit plugin adds `qiskit_regular_isometry`. `sparse_isometry_gf2x` is a deprecated alias, not a separate modern method. `identity_state_prep(num_qubits)` is a helper returning a circuit, not a registered algorithm that accepts a wavefunction.

The [common state-preparation interface](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/algorithms/state_preparation/state_preparation.py) takes a `Wavefunction` and returns a `Circuit`. That common signature must not hide output semantics: `alias_sampling` returns a coefficient-index PREPARE oracle with entangled workspace, not a clean arbitrary molecular trial state.

The common dense scattering helper obtains active determinants and their coefficients. The register width is `configuration_set.num_modes() * determinant.bits_per_mode()`. For the ordinary spin-half spatial-orbital configurations, there are two bits per spatial mode, alpha in the lower half and beta in the upper half. Synthetic compressed subspaces use one-bit-per-mode configurations. The helper places each coefficient at the integer formed from its determinant occupation bits; it does not assume coefficient-list position equals basis index. Duplicate indices are added in this helper.

That helper rejects empty coefficients, non-negligible imaginary parts, and widths above 32, and produces a real vector of length \(2^q\). The 32-qubit guard is a sanity boundary, not a feasible-memory recommendation: \(2^{32}\) eight-byte real amplitudes alone occupy 32 GiB. There is no cheap route through this helper for a generic sparse 56-qubit nitrogen state. QROM and alias sampling currently use the same dense scattering helper, so their names do not remove that classical allocation.

## Dense pure-state preparation

[`dense_pure_state`](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/algorithms/state_preparation/dense_pure_state.py) is the simplest native Q#-backed loader. It scatters a real wavefunction into the full occupation vector and constructs a `StatePreparationParams` record with a reversed `rowMap`, the amplitude vector, no expansion operations, and the register width. It returns a Q# callable and factory with `encoding="jordan-wigner"` and known `num_qubits`.

The repository’s [`StatePreparation.qs`](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/utils/qsharp/src/StatePreparation.qs) calls external `Std.StatePreparation.PreparePureStateD` on the subregister selected by `rowMap`. The reversal reconciles its vector/register convention with occupation-bit indices; it does not reverse the physical meaning of alpha and beta. Existing [dense-loader tests](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/tests/test_state_preparation_dense.py) compare prepared vectors, single determinants, and energies with sparse isometry.

The Python helper itself does not explicitly normalize the scattered vector despite broad wording in its method docstring. Loading behavior beyond the adapter is delegated to the standard-library operation. Supply a valid normalized, finite nonzero target rather than relying on a blanket claim of validation or repair. The native dense class adds no precision or transpiler settings of its own.

Conceptually, dense state loading can recursively split amplitude norm between branches. If a prefix has total probability \(p\), of which \(p_0\) lies in its zero-child, a Y rotation can allocate amplitudes proportional to \(\sqrt{p_0/p}\) and \(\sqrt{1-p_0/p}\). Repeating down a tree gives general amplitude loading. This explains the exponential generic information cost, but it is not a claim that we reviewed the external standard-library decomposition line by line.

An **isometry** preserves inner products from an input space into an output space. Preparing one state maps a known reference vector to that state; a quantum implementation extends that action to a unitary on a complete space, possibly with clean workspace. It does not need to implement an arbitrary unknown unitary supplied as an exponentially large matrix. Nevertheless a completely arbitrary dense target contains exponentially many independent parameters.

## The optional Qiskit regular isometry

[`qiskit_regular_isometry`](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/plugins/qiskit/regular_isometry.py) converts the wavefunction to a dense complex statevector, asks Qiskit’s `StatePreparation` to synthesize it, optionally transpiles it, and exports OpenQASM 3. This is an external Qiskit synthesis contribution, not a second native Q# implementation.

Its settings are `transpile=True`, `transpile_optimization_level=0`, and `basis_gates=["x","y","z","cx","cz","id","h","s","sdg","rz"]`. The adapter additionally applies repository passes `MergeZBasisRotations`, `RemoveZBasisOnZeroState`, and `SubstituteCliffordRz`. These are state-preparation-specific simplifications, not justification to remove such gates from arbitrary unknown-input unitaries.

The [conversion helper](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/plugins/qiskit/conversion.py) preserves complex coefficients and normalizes when its norm is above `1e-15`. It allocates \(2^q\) complex amplitudes without the native helper’s explicit 32-qubit guard. Duplicate determinant entries are assigned rather than summed on this path. Use unique, valid determinant support instead of assuming all loaders repair duplicates identically.

The circuit is returned with QASM and Jordan–Wigner encoding, but without an explicitly supplied `num_qubits` field in this adapter. Conversion/building consumers must handle that declared-width limitation rather than assume every method returns identical metadata. Dense complex support here does not imply sparse isometry accepts complex input: the sparse parent checks for real coefficients before delegating.

## Sparse isometry: compress the support, load, then expand

The default [`sparse_isometry`](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/algorithms/state_preparation/sparse_isometry.py) exploits the bit patterns of the supported determinants. Build a binary matrix with one row per qubit and one column per determinant. The entries say which occupation bits are one. Coefficients remain aligned with columns throughout; ordinary floating-point row reduction is not the transformation being applied to amplitudes.

A CNOT from row \(u\) to row \(v\) replaces row \(v\) by its XOR with row \(u\). An X gate complements one row. Both are reversible operations on bitstrings. If two rows have identical patterns, a CNOT zeros one. If a row is all ones across the support, X zeros it. Binary Gaussian elimination then finds independent rows and records the CNOTs that eliminated dependent patterns. Removing zero rows at this stage means those qubits are known to be zero **in the transformed supported state**, not that their physical modes have been deleted from the Hamiltonian.

The implementation performs duplicate-row removal, all-ones-row removal, GF(2) elimination, zero-row removal, and an additional guarded diagonal/staircase reduction when its structural test permits it. It tracks the retained row map and operations. This is affine compression: the X translations augment linear XOR transformations. The retained width is related to the affine span of the supported bitstrings, not simply the real matrix rank or the number of determinants.

For a hand-worked example, suppose support contains only the increasing-qubit-order bit vectors `0000` and `1111`, with amplitudes \(\sqrt{0.8}\) and \(\sqrt{0.2}\). All four rows of the support matrix agree. Three CNOT eliminations can leave one independent row. Load \(\sqrt{0.8}|0\rangle+\sqrt{0.2}|1\rangle\) on that qubit and reverse the three CNOTs to reconstruct the four-qubit state. The final circuit still acts on four problem qubits. The saving was a one-qubit dense-loading stage, not an alteration of the chemical model.

For \(d\) determinants and retained width \(r\), the support matrix costs roughly \(qd\) bits in principle, though the Python implementation uses byte-valued arrays and other bookkeeping. The reduced amplitude vector has \(2^r\) entries. Even modest support can have a large affine span, so this method is not guaranteed cheap merely because \(d\ll2^q\).

The ordinary path places each coefficient at its reduced determinant index, normalizes the nonzero reduced vector, constructs a synthetic one-bit-per-mode wavefunction, and delegates dense loading. Its `dense_state_prep` setting defaults to `AlgorithmRef("state_prep","dense_pure_state")`. It then reverses the recorded GF(2)+X operations to restore the original register. Default expansion gates are Clifford X and CNOT gates. Creating the synthetic wavefunction is a software representation of compressed bits, not inventing a new molecule or a smaller physical active space.

The native nested dense factory is composed in Q# using its parameters; this avoids a callable-capturing-callable limitation in adaptive code generation. A nested Qiskit loader is composed via QASM/QIR and Qiskit. A custom child circuit carrying only a prebuilt Q# circuit object is rejected if it cannot enter either composition path. The parent’s [`dense_state_prep` contract](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/algorithms/state_preparation/sparse_isometry.py) is therefore more specific than “any algorithm in the family works.”

There is no exposed support-pruning threshold in the sparse parent. Its comments mentioning filtering do not implement automatic top-determinant selection. Remove zero coefficients, merge duplicate determinants, select the desired support, and verify normalization beforehand. Duplicate columns can otherwise overwrite coefficients during reduced scattering; an all-identical multi-entry support can produce a zero-rank error. Non-negligible imaginary coefficients are rejected even when the nested loader is complex-capable.

For one determinant, the sparse method bypasses dense loading and applies X on each occupied bit. The actual Q# body uses `bitStrings[i]` on `qs[i]`; a contradictory helper docstring about reversing the input does not describe its behavior. A normalized single determinant’s overall sign is physically global and is not retained by this shortcut. A zero-amplitude “single determinant” should not be supplied as a state. Our six-electron occupation determinant needs six X gates on twelve problem qubits, not six qubits or a zero-cost identity operation.

## Optional binary compression and its non-Clifford price

The sparse method’s `binary_encoding=False` setting enables a second, best-effort route when changed to true. It first obtains forward row-echelon form, then compares the remaining row count with \(\lceil\log_2 d\rceil\). It takes this route only when the latter is strictly smaller; otherwise it falls back to ordinary sparse isometry.

This compression is **not** Bravyi–Kitaev fermionic encoding. It assigns compact binary labels to the finite supported determinant set, loads amplitudes on those labels, and expands back to the same Jordan–Wigner output register.

The [_binary_encoding_utils implementation](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/algorithms/state_preparation/_binary_encoding_utils.py) first compresses pivot columns from a unary/staircase form to binary. It then processes remaining columns in batches, using Toffoli operations and lookup blocks called PUI blocks in the code to clear sparse indicator rows. A **Toffoli** is a NOT controlled on two bits; unlike the default X/CNOT expansion it is not Clifford. A recorded bijection keeps compact labels matched to original coefficient columns. The synthesizer verifies that non-dense rows are zero and its mapping covers the columns.

`include_negative_controls=True` permits controls conditioned on a bit being zero as well as one. The implementation promotes fixed control patterns and can split batches into power-of-two chunks when its lookup-cost estimate favors that. Smaller dense width can therefore trade arbitrary rotations for Toffolis, lookup tables, and workspace. The winner depends on those costs, not only on a new label count.

`measurement_based_uncompute=False` is the other control. When enabled, selected lookup branches use `Std.Intrinsic.AND` and its adjoint instead of a fully unitary Toffoli uncompute. The [Q# binary expansion](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/utils/qsharp/src/BinaryEncoding.qs) borrows known-zero qubits from an ancilla pool formed from unused transformed rows, returns them clean, and allocates helpers if that pool is insufficient. Mid-circuit measurement with corrective feedback can reduce non-Clifford work while preserving the intended logical action. It is not measurement of the final chemical occupations.

The target profile matters. The adapter warns that a Base profile cannot use the intended measurement-based path and relies on unitary uncomputation there; adaptive profiles can lower it to measurement and feedforward. Existing [binary-state tests](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/tests/test_state_preparation_binary_encoding.py) check reduced CCZ counts, added measurements, QIR/Qiskit control flow, and output-state independence from measurement outcome. Those are source checks, not new measured savings here. The binary composition advertises `Adj`, not automatically `Ctl`; adjointable does not mean every arbitrary outer-controlled use is supported.

`create_dense(wavefunction)` returns only the dense-loading component embedded in the same full-width register, omitting expansion. It is for component costing. It is generally **not** the full requested trial state. Subtracting additive logical gate counts from matched full/component estimates can isolate expansion work; subtracting total physical resource estimates need not have the same meaning because scheduling and error correction are nonlinear.

## QROM: loading a fixed classical table into rotations

**QROM**, quantum read-only memory in these circuits, means a reversible circuit implementing a lookup table known at circuit construction. It is not a free, universal hardware QRAM device. A lookup correlates an address register with table output, and that temporary output must subsequently be uncomputed.

[`state_prep/qrom`](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/algorithms/state_preparation/qrom_state_prep.py) is a real signed-amplitude loader. It scatters the target over the full \(2^q\) occupation vector, checks finite coefficients and squared coefficients and rejects an all-zero vector, then uses `rotation_bit_precision=10` by default, with the exposed range 1–30. Rescale extreme data rather than regard these checks as a complete safeguard against every possible overflow in a later probability sum.

The [Q# implementation](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/utils/qsharp/src/QROMStatePrep.qs) explicitly implements a Shende–Bullock–Markov-style probability tree. At a node with total mass \(p\) and left-child mass \(p_L\), it forms \(\theta=\arccos\sqrt{p_L/p}\); zero-mass nodes get zero angle. It needs a rotation \(R_y(2\theta)\) to distribute the amplitudes. For precision \(b\), it stores the integer

\[
x=\operatorname{Round}\!\left(\frac{2^b\theta}{2\pi}\right)\bmod2^b.
\]

The phase-gradient rotation routine implements \(R_y(4\pi x/2^b)\), matching that convention. Thus quantization changes the prepared amplitudes by a controlled numerical approximation, not by an electron-space reduction. Ignoring wrap-equivalent angles, nearest rounding bounds the individual rotation-angle error by \(2\pi/2^b\); errors across many rotations still require a total state/error budget.

The circuit fixes the most significant occupation bit first, using already fixed prefix bits as lookup addresses for later rotations. It reverses slices explicitly where the little-endian table-lookup interface requires it. Magnitude rotations alone would lose coefficient signs, so a separate sign table is loaded into a scratch sign bit; Z applies the required sign phase, and the lookup is undone. The [signed-amplitude tests](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/tests/test_state_preparation_qrom.py) check mixed signs, all-negative vectors, zero entries, and determinant-index placement.

The repository’s [`SelectSwap.qs`](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/utils/qsharp/src/SelectSwap.qs) splits an address between selection and a swap network. More swap bits mean more parallel scratch data words but fewer select addresses. Passing −1 asks the helper to search its internal integer cost model for a split. The no-swap select model is linear in table length; swap, scratch, and erasure terms change the tradeoff. Do not turn an abstract lookup into a unit-cost oracle, or interpret that internal logical estimate as a measured hardware runtime. The total fixed table and the bare Python dense allocation remain potentially exponential.

## Reusable phase-gradient ancillas

A **phase-gradient state** of \(b\) qubits has amplitudes

\[
|\phi_b\rangle=2^{-b/2}\sum_{k=0}^{2^b-1}e^{-2\pi i k/2^b}|k\rangle.
\]

Adding integer \(x\) modulo \(2^b\) leaves this state unchanged except for phase \(e^{2\pi i x/2^b}\). Condition the addition appropriately and that phase is kicked back onto a target qubit. This makes an integer-valued angle register useful for rotations without synthesizing a completely new arbitrary rotation at every lookup.

[`PhaseGradient.qs`](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/utils/qsharp/src/PhaseGradient.qs) prepares the gradient using X and an adjoint QFT, performs addition through an external QDK arithmetic primitive, and uses basis changes for its Y rotation. QFT means quantum Fourier transform. The standard-library QFT/adder internals are external-team contributions, not code we claim to have reviewed here. The routine’s documented cost model excludes preparing and unpreparing the gradient itself.

The QROM setting `allocate_phase_gradient=True` makes the callable self-contained: allocate the gradient internally, prepare it, use it, and unprepare it. Its declared `num_qubits` is then just the target-state width, and `num_phase_gradient_ancillas` is zero because no gradient is expected from the caller. Internal allocation still contributes to actual peak logical width.

With `allocate_phase_gradient=False`, the callable expects `[state register | gradient register]`. It declares `num_qubits=q+b` and `CircuitMetadata(num_phase_gradient_ancillas=b)`. The caller must initialize the trailing gradient and preserve it across compatible subroutines. It must not include that nonzero gradient in a reflection intended to test whether ordinary ancillas are all zero. Metadata identifies the sharing contract; it does not prepare the gradient by itself.

There is a subtle export boundary: the associated standalone factory always allocates and prepares its own gradient because an exported circuit has no caller to lend one. Costing the factory and costing a composed callable with a shared gradient answer different questions. Existing tests explicitly distinguish those declared widths. Reuse can amortize preparation **within one composed execution**, not across independently initialized jobs without a maintained quantum state.

Precision limits such as 30 bits are engineering bounds. A comment that an angle grid is “below chemical accuracy” is not a proof about total molecular energy error, long controlled evolutions, or noisy preparation.

## Alias sampling: the same marginal distribution is not the same pure state

[`state_prep/alias_sampling`](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/algorithms/state_preparation/alias_sampling.py) accepts real **nonnegative** wavefunction coefficients, squares them after dense determinant-index scattering, and passes those weights to Q#. It rejects negatives, non-finite values, overflow on squaring/summing, and all-zero weights. Its `bits_precision` defaults to 10 with range 1–30.

For normalized target amplitudes \(c_j\ge0\), the target marginal probabilities are \(p_j=c_j^2\). The result is generally

\[
\sum_j\sqrt{\widetilde p_j}|j\rangle|g_j\rangle,
\]

where \(g_j\) is workspace or **garbage** correlated with the index, and \(\widetilde p\) approximates \(p\). Tracing out that workspace does not in general yield the pure state \(\sum_j c_j|j\rangle\). Correct occupation probabilities alone cannot establish correct interference or ground-state overlap.

The [Q# alias implementation](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/utils/qsharp/src/AliasSamplingStatePrep.qs) first discretizes probabilities into \(L2^\mu\) integer units, where \(L\) is table length and \(\mu\) the precision. It floors each allocation, distributes leftover units by largest fractional remainder with deterministic ties, and ensures zero weights do not gain population. The resulting per-entry probability error is bounded by \(1/(L2^\mu)\) in the intended discretization; the [tests](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/tests/test_state_preparation_alias.py) check that bound and zero-weight behavior.

Walker-style alias construction then assigns each table position a keep count and an alternate index, pairing underfilled and overfilled probability bars. Coherently, the circuit prepares a uniform index, a uniform \(\mu\)-bit comparison value, loads `(keep, alternate)` through QROM, compares the value with keep, and conditionally swaps the index with its alternate. The flag/comparison/lookup registers retain the information needed for reversibility. They are not harmless classical notes one can discard.

Named registers total \(2\lceil\log_2L\rceil+2\mu+1\): index, alternate storage, comparison bits, keep storage, and flag. Lookup/swap implementations can add more scratch. In this Python wrapper, \(L\) is the full scattered occupation-vector length, not automatically the sparse determinant count. A sparse two-entry wavefunction at indices 0 and 3 is explicitly scattered to four weights in a source test.

This is useful as **PREPARE** in a linear-combination-of-unitaries, or LCU, block-encoding construction, where SELECT acts on the index and PREPARE-adjoint is part of the coherent structure. Even there, unpreparing the garbage is not a claim that every branch deterministically returns all workspace to zero after an arbitrary SELECT. The useful zero-ancilla block is extracted by the construction we will derive later.

To load an LCU distribution proportional to \(|w_j|\) through this particular Python wavefunction adapter, the supplied nonnegative amplitudes must be proportional to \(\sqrt{|w_j|}\), because it squares them. Signs of LCU coefficients belong in the matching SELECT action. Separate conditional Q# alias helpers do carry sign bits and conditional/free-rider data for specialized builders; that does not lift the public `alias_sampling` adapter’s nonnegative restriction.

The Q# code uses a phase-safe ancilla-free uniform preparation for non-power-of-two lists, because PREPARE-adjoint must work after phases have been imprinted, not only on the pristine state it first prepared. This is another reason an apparently cheaper preparation shortcut may be invalid as a reusable coherent subroutine.

## Identity preparation, useful code patterns, and complete cost

[`identity_state_prep(q)`](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/algorithms/state_preparation/identity.py) wraps an empty repeated Pauli sequence. It leaves the input register unchanged. On fresh zeroed qubits it therefore leaves the vacuum bitstring, not a six-electron nitrogen determinant. It does not inspect a wavefunction or load a desired state for free. The helper currently supplies neither encoding metadata nor an explicit `num_qubits` field to `Circuit`, although its factory records target indices. Builders requiring those fields cannot infer a chemistry contract from the word identity.

For the continuing untapered Jordan–Wigner story, a source-backed, unexecuted pattern is:

```python
from qdk_chemistry.algorithms import create
from qdk_chemistry.data import AlgorithmRef, Ansatz

# trial_wavefunction contains valid real CI coefficients in the selected orbitals.
trial = Ansatz(active_hamiltonian, trial_wavefunction)
prep = create(
    "state_prep", "sparse_isometry",
    dense_state_prep=AlgorithmRef("state_prep", "dense_pure_state"),
    binary_encoding=False,
)
state_circuit = prep.run(trial.get_wavefunction())
```

Changing the nested reference to `qiskit_regular_isometry` requires the optional plugin; its transpilation knobs belong on that nested reference, not a fictional global state-preparation setting. Enabling binary compression is a tradeoff to evaluate on actual support, not an unconditional improvement. A QROM alternative is `create("state_prep", "qrom", rotation_bit_precision=10).run(trial_wavefunction)` under its real-amplitude and dense-width limits.

A fair comparison includes classical CI/support selection, preprocessing tables and bit matrices, dense loading, expansion, sign/phase handling, all scratch allocation, gradient initialization, adjoints, and every repetition the enclosing algorithm requires. It also distinguishes logical rotations from the fault-tolerant gates used to approximate them. Adding determinants may improve overlap enough to save trials, or may cost more in loading than it saves. No single overlap or gate-count number establishes that balance.

For our nitrogen baseline the state register remains twelve qubits even when a sparse loader uses a much smaller dense subregister. A tapered ten-qubit operator from Lesson 15 cannot consume this state circuit without a matched transformation. A localizer’s one-determinant carrier cannot replace the correlated CI state by retaining its old filename. And a coefficient-loading oracle with entangled garbage cannot stand in for a clean trial state because its marginal histogram looks right.

We can now ask the next engineering question precisely: given a state-preparation subroutine and a mathematical Hamiltonian operation in compatible representations, how does the repository assemble them into actual circuits and interchange formats? That is the subject of Lesson 17.