Exact numerical simulation

analysis/README.md

@@ -284,7 +284,79 @@ There are many details of the algorithm that may be worth analyzing. The availab
     - `"largest"` gives [10, 5, ...] (most positive)
     - This is NOT based on magnitude (which would give [0, 5, -5, ...])
 
-#### Numerical Simulation
+#### Error Analysis
+
+- **Error Analysis**: Enables three types of error metrics comparing exact and approximate representations. Each error type is independently enabled by setting its corresponding configuration parameter.
+
+  **Configuration parameters**:
+
+  - **`enable_eigenvalue_errors`**: Enable eigenvalue error comparison (default: False, disabled)
+    - Set to `True` to compute errors for ALL eigenvalues from the eigendecomposition
+    - The number of eigenvalues compared is determined by `num_eigenvalues` setting
+    - Compares all eigenvalues computed in both exact and approximate eigendecompositions
+    - Requires `eigendecomposition_matrices = "both"` to ensure both decompositions are computed
+    - **Best for**: Validating eigenvalue accuracy across the spectrum
+
+  - **`error_matrix_norms`**: Which matrix norms to compute (default: None, disabled)
+    - Single string: `"frobenius"` or `"spectral"`
+    - List: `["frobenius", "spectral"]` for both
+    - **Frobenius norm**: Element-wise difference, fast to compute
+      - ||H_exact - H_approx||_F = sqrt(sum of squares of all elements)
+      - Good for quick comparisons
+    - **Spectral norm**: Worst-case effect on any quantum state, physically meaningful
+      - ||H_exact - H_approx||_2 = largest singular value
+      - More expensive to compute, especially for large systems
+    - **For large systems (>15 qubits)**: Uses matrix-free computation
+      - Frobenius: Requires 2^N matrix-vector products (minutes for 20 qubits)
+      - Spectral: Uses power iteration (can take longer)
+      - Progress warnings displayed during computation
+    - **Best for**: Physical bounds on algorithm error
+
+  - **`error_state_inputs`**: State files for state-dependent errors (default: None, disabled)
+    - Single filename (string): `"ground_state.npy"`
+    - Multiple filenames (list): `["state1.npy", "state2.npy"]`
+    - Compares: ||H_exact|ψ⟩ - H_approx|ψ⟩||
+    - **Best-scaling error metric** for large systems
+      - Only requires O(2^N) memory (state vectors), not O(2^(2N)) (matrices)
+      - Can reach 30 qubits
+      - Fast: just applies operators to states
+    - **Best for**: Error on specific physically relevant states
+
+  **System size guidance**:
+  - **≤15 qubits**: All error types fast (dense matrix operations)
+  - **16-20 qubits**: Matrix norms slow (matrix-free), state errors still fast
+  - **20-30 qubits**: Use eigenvalue + state errors; avoid matrix norms unless necessary
+  - **Production recommendation**: Eigenvalue errors (k=1-5) + state errors for best performance
+
+  **Output file**: `error_analysis.npz` containing all computed error metrics
+
+  **Examples**:
+  ```python
+  # Eigenvalue error comparison (compares all eigenvalues from eigendecomposition)
+  analysis.enable_eigenvalue_errors = True
+
+  # Matrix norm errors (physical bounds)
+  analysis.error_matrix_norms = "frobenius"  # Fast
+  # or
+  analysis.error_matrix_norms = ["frobenius", "spectral"]  # Both
+
+  # State-dependent errors (best scaling)
+  analysis.error_state_inputs = "ground_state.npy"
+  # or multiple states
+  analysis.error_state_inputs = ["ground.npy", "excited.npy"]
+
+  # All three error types together (comprehensive validation)
+  analysis.enable_eigenvalue_errors = True
+  analysis.error_matrix_norms = "frobenius"
+  analysis.error_state_inputs = ["ground.npy", "excited.npy"]
+  ```
+
+  **When to use each error type**:
+  - **Eigenvalue errors**: Use when you want to validate all eigenvalues computed in the eigendecomposition
+  - **Matrix norm errors**: Use for physical bounds, but expensive for large systems
+  - **State errors**: Use for specific physically relevant states, scales best
+
+#### Numerical Simulation (Approximate)
 
 - **Numerical Simulation**: Setting `analysis.numerical_simulation_inputs` to one or more state
   vector files will apply the algorithm to the input state(s) via numerical simulation, producing
@@ -321,7 +393,56 @@ There are many details of the algorithm that may be worth analyzing. The availab
   psi[0] = 1.0
   np.save("initial_state.npy", psi)
   ```
+
+#### Exact Numerical Simulation
+
+- **Exact Numerical Simulation**: Setting `analysis.exact_simulation_inputs` to one or more state
+  vector files will apply the **exact Hamiltonian matrix** (without Trotter or other approximations)
+  to the input state(s), producing output state(s).
+
+  **Purpose**: Compare exact time evolution with approximate algorithm results
+
+  **Input format**: Same as approximate simulation (NumPy `.npy` format). Input can be:
+  - Single filename (string): `"initial_state.npy"`
+  - Multiple filenames (list): `["state1.npy", "state2.npy", "state3.npy"]`
+
+  **Output naming**: Automatic suffix `_exact_final` is added to input filename:
+  - `initial_state.npy` → `initial_state_exact_final.npy`
+  - This distinguishes from approximate simulation output (`_final`)
+
+  **Scaling**: Can reach 25-30 qubits using matrix-free operators from Branch 1
+  - For large systems, uses `PauliStringOperator` (matrix-free)
+  - Memory: O(2^N) for state vectors, not O(2^(2N)) for matrices
+  - Computation time depends on number of Pauli terms in Hamiltonian
+
+  **Example**:
+  ```python
+  # Single state exact simulation
+  analysis.exact_simulation_inputs = "initial_state.npy"
+
+  # Multiple states
+  analysis.exact_simulation_inputs = [
+      "ground_state.npy",
+      "excited_state.npy",
+      "superposition.npy"
+  ]
+  ```
+
+  **Comparing exact vs approximate**:
+  You can run both simulations on the same input states:
+  ```python
+  # Run both approximate and exact simulations
+  analysis.numerical_simulation_inputs = "initial_state.npy"
+  analysis.exact_simulation_inputs = "initial_state.npy"
 
+  # Produces:
+  #   initial_state_final.npy       (approximate result)
+  #   initial_state_exact_final.npy (exact result)
+  # Then manually compare the output states
+  ```
+
+  **Note**: State vectors must be compatible with Hamiltonian dimension (2^N)
+
 ## Generated Files
 
 The script will print a log both to the screen and to a logfile. It also generates output files
@@ -332,8 +453,14 @@ depending on the analyses requested:
   results. Shows final results but not intermediate values.
 - **Matrix file** (if `matrix_output_file` specified): Unitary matrix in specified format (`.npz`,
   `.h5`, or `.txt`)
+- **Exact matrix file** (if `exact_matrix_output_file` specified): Exact Hamiltonian matrix
+- **Eigendecomposition files** (if `num_eigenvalues` > 0): `exact_eigendecomposition.npz` and/or
+  `approximate_eigendecomposition.npz`
+- **Error analysis file** (if any error analysis enabled): `error_analysis.npz`
 - **Final state files** (if `numerical_simulation_inputs` specified): Evolved quantum states with
   `_final` suffix (e.g., `initial_state.npy` → `initial_state_final.npy`)
+- **Exact final state files** (if `exact_simulation_inputs` specified): Exactly evolved states with
+  `_exact_final` suffix (e.g., `initial_state.npy` → `initial_state_exact_final.npy`)
 
 The logfile is typically most useful for understanding the analysis process and intermediate values.