Add additional analyses for studying Trotter error

[BrendanKKrueger]

As a first pass (because I'm about to leave on vacation), I put together a plan with Claude and then let Claude run.

Merge requests:

• Comprehensive configuration example #26
• Compute the exact matrix #27
• Perform eigendecomposition #28
• Compute error metrics #29
• Exact numerical simulation #30

Initial Prompt

Please develop a plan to add the following capabilities to the analysis tool:

• Compute the eigendecomposition of the unitary matrix equivalent to the algorithm.
• Compute the exact matrix form for the algorithm (that is, no Trotter or double-factorization, but the exact version).
• Compute the eigendecomposition of the exact matrix.
• Perform numerical simulation using the exact matrix.
• Compute the eigenvalue error: the difference between the eigenvalue from the exact matrix and the eigenvalue from the unitary matrix equivalent to the algorithm. This will probably compare the lowest eigenvalues.
• Compute the matrix-norm error: the norm of the difference between the unitary equivalent of the algorithm and the exact matrix.
• Compute the state-dependent error: compare the final state from numerical simulation using the unitary equivalent to the algorithm and the final state from numerical simulation using the exact matrix.

Use the branch bkk_numerical as a starting point on which to build new changes.

The analyses discussed here have significant overlap in the things they need to compute. As with the bkk_numerical branch, any shared results (especially intermediate results from expensive calculations) should be computed once and then re-used as appropriate. Where appropriate, memory should be cleaned up once intermediate results that use significant memory are no longer needed.

These tasks should be appropriately distributed across multiple git merge requests for "bite size" reviews. Branches should either be independent or should build on each other as appropriate, but any dependencies between branches/merge requests needs to be clearly indicated.

All appropriate testing should be included (unit and integration). All appropriate documentation (README files, code comments, etc) updates should be made. Intermediate documentation explaining your process should not be added to the repository. Examples should be expanded and/or new examples should be added -- see analysis/examples/config2.py for an example of an example configuration file that is intended to exercise all analyses. The config2.py file may need to be updated to reflect all current analyses and also renamed before beginning the above tasks.

All tests (unit and integration) should be added to both the GitHub and GitLab CI pipelines. All tests should be run for each deliverable (that is, each branch that will be used for a merge request).

I will likely be unable to work through all these tasks before leaving on travel, so save the plan and the current status (at all reasonable break points, such as completing a branch, etc) to a file so that it can be picked up and continued later.

Claude's Plan

Context

The qhat analysis tool currently supports:

• Computing unitary matrices from quantum algorithms (with Trotter/double-factorization approximations)
• Performing numerical simulations using these unitary matrices
• Resource estimation

This plan adds capabilities to compute exact Hamiltonian matrices (without Trotter or approximation), their eigendecompositions, and various error metrics comparing the approximate algorithm to the exact solution. These features are critical for validating algorithm accuracy and understanding the impact of approximations.

The work builds on the bkk_numerical branch, which introduced:

• Unitary matrix computation via _compute_unitary_matrix()
• Numerical simulation with numerical_simulation()
• Shared intermediate result pattern (compute once, reuse)
• Analysis framework in analysis/analysis.py

Recent main branch updates (merged after bkk_numerical diverged):

• Major logging refactor: Centralized logging architecture with common/logging_utils.py
• All modules now use standard logging.getLogger(__name__) pattern
• Removed mock_config.py from tests in favor of conftest.py with logging fixtures
• config_types.py: Removed inline logging configuration (now handled by logging_utils)
• analysis.py: All functions now use module-level logger instead of config_general.log_*()
• These changes will need to be merged into bkk_numerical before proceeding with implementation

Architecture Overview

Current Structure

• analysis/analysis.py: Main analysis dispatcher with analyze_algorithm()
• analysis/config_types.py: Configuration classes (AnalysisConfiguration, etc.)
• analysis/hamiltonian.py: Hamiltonian class with Pauli string generation
• analysis/unitary.py: Unitary encoding methods (Trotter, double-factorization)
• analysis/file_io.py: Matrix and state I/O utilities
• Shared results pattern: Expensive computations (like unitary matrix) computed once and passed to multiple analyses

Required Additions

1. Exact matrix computation: Convert Hamiltonian Pauli strings directly to matrix form
2. Eigendecomposition: Compute eigenvalues/eigenvectors for both exact and unitary matrices
3. Error metrics: Three types of errors comparing exact vs approximate
4. Configuration: Extend AnalysisConfiguration to request new analyses
5. Testing: Unit tests for matrix operations, integration tests for full pipeline
6. CI/CD: Add new tests to GitHub Actions and GitLab CI

Implementation Plan: Branch Structure

IMPORTANT: Merge Strategy

Before starting Branch 0, the bkk_numerical branch must be updated with the latest changes from main. The main branch has undergone a significant logging refactor that affects:

• How modules obtain loggers (now logging.getLogger(__name__))
• Test infrastructure (removed mock_config.py, added conftest.py)
• Configuration classes (removed inline logging setup)

Recommended approach:

1. Merge main into bkk_numerical: git merge main
2. Resolve conflicts in files that both branches modified:
   • analysis/analysis.py: bkk_numerical added numerical simulation; main refactored logging
   • analysis/config_types.py: bkk_numerical added numerical_simulation_inputs; main refactored logging
   • analysis/file_io.py: bkk_numerical added state I/O; main refactored logging
   • analysis/hamiltonian.py: Both branches have minor changes
   • analysis/tests/: bkk_numerical added new tests; main added conftest.py and removed mock_config.py
3. Update bkk_numerical code to use new logging patterns:
   • Replace config_general.log() with logger.info()
   • Replace config_general.log_verbose() with logger.verbose()
   • Replace config_general.log_debug() with logger.debug()
   • Update test files to use conftest.py instead of mock_config.py
4. Run all tests to verify compatibility: python3.11 -m pytest analysis/tests/ -v
5. Then proceed with Branch 0

Expected conflicts:

• All logging calls in analysis.py functions (numerical_simulation, output_unitary_matrix)
• Test imports in test_numerical_simulation.py (references mock_config.py which was removed)
• Minor formatting differences in other files

Files unique to bkk_numerical (no conflicts, just keep them):

• analysis/tests/test_numerical_simulation.py
• analysis/tests/test_analysis_integration.py
• Enhanced analysis/README.md sections

Alternative approach (if merge is complex):

• Rebase bkk_numerical onto current main: git rebase main
• This gives a cleaner history but requires more careful conflict resolution
• Same logging pattern updates will be needed

The work is organized into 5 sequential branches (each builds on the updated bkk_numerical):

Branch 0: bkk_config_full_analysis (Baseline Configuration)

Purpose: Establish comprehensive example configuration that exercises all current QHAT analyses

Builds on: bkk_numerical

Files to modify:

• Rename analysis/examples/config2.py → analysis/examples/config_full_analysis.py
• Update the file to include ALL currently available analyses:
  • Resource estimation (pyLIQTR)
  • Matrix output (unitary matrix from algorithm)
  • Numerical simulation
• Update analysis/README.md to reference config_full_analysis.py instead of config2.py
• Update any other documentation that references config2.py

Key implementation details:

• This branch adds NO new features - just improves the example
• The configuration should demonstrate:
  • All resource estimation options
  • Matrix output to different formats (npz/h5/txt)
  • Numerical simulation with multiple input states
  • Comments explaining what each analysis does and when to use it
• Acts as baseline: subsequent branches will add new analyses to this example

Testing:

• Run config_full_analysis.py and verify all existing analyses work
• Verify output files are created as expected
• No new tests needed (existing tests cover functionality)

Merge criteria:

• Example runs successfully on test system
• README updated
• All references to config2.py updated
• CI/CD passes (no functional changes, just rename/reorganization)

---

Branch 1: bkk_exact_matrix (Foundation)

Purpose: Compute exact Hamiltonian matrices without approximations

Builds on: bkk_config_full_analysis

Files to modify:

• analysis/analysis.py: Add _compute_exact_matrix() helper and exact_matrix_output() analysis
• analysis/config_types.py: Add exact_matrix_output_file field to AnalysisConfiguration
• analysis/hamiltonian.py: Add to_matrix() method to Hamiltonian class
• analysis/README.md: Document new exact matrix capability
• analysis/examples/config_full_analysis.py: Add example of exact matrix output
• .github/workflows/python-app.yml: Ensure tests run
• .gitlab-ci.yml: Ensure tests run

New files:

• analysis/matrix_operations.py: Utility functions for Pauli string operations
  • pauli_char_to_matrix(): Convert 'I', 'X', 'Y', 'Z' to 2x2 matrices (for small systems)
  • pauli_string_to_matrix(): Build full matrix from string (for small systems)
  • pauli_dict_to_matrix(): Sum weighted Pauli strings into dense matrix (for small systems)
  • PauliStringOperator class: Matrix-free linear operator for large systems
    • Stores Hamiltonian as Pauli strings (sparse representation)
    • Implements matvec(v) to apply operator to vector without materializing matrix
    • Compatible with scipy sparse eigensolvers
    • Scales to 25-30+ qubits (limited by state vector size, not matrix size)
• analysis/tests/test_matrix_operations.py: Unit tests for matrix utilities
• analysis/tests/test_exact_matrix.py: Integration tests for exact matrix computation

Key implementation details:

• Reuse hamiltonian.get_all_pauli_strings(return_as='strings') to get coefficient dict

• Two-path design based on system size:

  Small systems (≤15 qubits): Dense path

  • Materialize full dense matrix using existing pauli_dict_to_matrix()
  • Store as numpy array
  • Fast, simple, reusable

  Large systems (>15 qubits): Matrix-free path

  • Create PauliStringOperator that stores Pauli strings
  • Implements matvec(v) for matrix-vector products without materializing matrix
  • Compatible with scipy LinearOperator interface
  • Memory: O(num_pauli_terms + 2^N) instead of O(2^(2N))
  • Can reach 25-30+ qubits

• Automatic selection: Check dimension, use dense if ≤32768 (15 qubits), else matrix-free

• User override: Optional force_dense or force_sparse parameters

• Matrix output to file:

  • Dense: save directly via save_matrix()
  • Matrix-free: Either refuse (dimension too large) or save Pauli string representation

• Key insight: Different analyses may need different representations - some cache, some compute on-the-fly

Testing:

• Unit: Test Pauli conversion functions with known matrices
• Unit: Test matrix properties (Hermiticity, eigenvalue bounds match 1-norm)
• Integration: Test full pipeline with small Hamiltonian (2-3 qubits)
• Integration: Verify exact matrix matches expected form for test cases

Merge criteria:

• All unit and integration tests pass
• GitHub Actions and GitLab CI green
• README updated with exact matrix documentation
• Example config demonstrates usage

---

Branch 2: bkk_eigendecomposition (Eigenanalysis)

Purpose: Compute eigendecompositions of exact and unitary matrices

Builds on: bkk_exact_matrix

Files to modify:

• analysis/analysis.py: Add eigendecomposition_analysis() function
• analysis/config_types.py: Add eigendecomposition config fields to AnalysisConfiguration:
  • num_eigenvalues: int or str (number of eigenvalues to compute, default=0)
    • 0 = disabled (no eigendecomposition)
    • positive int = compute that many eigenvalues (e.g., 5 for ground state + 4 excited states)
    • "all" = compute full eigendecomposition (all eigenvalues and eigenvectors)
  • eigendecomposition_matrices: str ('exact', 'approximate', or 'both', default='approximate')
    • 'approximate' = only eigendecompose the algorithm's unitary (default - avoids extra computation)
    • 'exact' = only eigendecompose the exact Hamiltonian matrix
    • 'both' = eigendecompose both (useful for comparison, but computes both matrices)
  • which_eigenvalues: str ('smallest', 'largest', or 'both', default='smallest')
    • 'smallest' = algebraically smallest (most negative, ground state for Hamiltonians)
    • 'largest' = algebraically largest (most positive)
    • NOT based on magnitude/absolute value
    • Ignored when num_eigenvalues = "all" (full decomposition includes all)
• analysis/file_io.py: Add save_eigendecomposition() and load_eigendecomposition()
• analysis/README.md: Document eigendecomposition analysis
• analysis/examples/config_full_analysis.py: Add eigendecomposition configuration

New files:

• analysis/tests/test_eigendecomposition.py: Tests for eigendecomposition functionality

Key implementation details:

• Eigendecomposition triggered when num_eigenvalues != 0 (either positive int or "all")
• User controls which matrices to compute/decompose via eigendecomposition_matrices:
  • 'approximate' (default): Only compute approximate matrix and its eigendecomposition → approximate_eigendecomposition.npz
  • 'exact': Only compute exact matrix and its eigendecomposition → exact_eigendecomposition.npz
  • 'both': Compute both matrices and both eigendecompositions → two output files

Matrix representation strategy:

• For eigendecomposition, we don't always need to materialize the full matrix
• scipy.sparse.linalg.eigsh() accepts either:
  • A matrix (dense or sparse)
  • A LinearOperator with matvec(v) method
• Use matrix-free operator for large systems (dimension > 32768)
  • Pass PauliStringOperator or approximate algorithm's operator directly to eigsh
  • Computes eigenvalues without ever materializing the full matrix
  • Enables 25-30 qubit eigendecomposition

Eigendecomposition modes:

• Full eigendecomposition (num_eigenvalues = "all"):

  • Use scipy.linalg.eigh() to compute all eigenvalues/eigenvectors
  • Requires dense matrix - automatically refuse if dimension > 32768
  • Suitable for systems ≤15 qubits only
  • which_eigenvalues parameter ignored (all eigenvalues computed anyway)
  • Clear error message: "Full eigendecomposition not supported for >15 qubits, use num_eigenvalues=k"

• Partial eigendecomposition (num_eigenvalues > 0):

  • Use sparse methods: scipy.sparse.linalg.eigsh() with k parameter
  • Works with matrix-free operators - scales to 25-30+ qubits
  • Respects which_eigenvalues setting
  • Can compute k << dimension efficiently even for very large matrices
  • This is the primary production path

• Implementation: Check if isinstance(num_eigenvalues, str) and num_eigenvalues.lower() == "all"

• Memory management: eigsh only stores k eigenvectors (not full basis)

• Output format for each npz file:

  • eigenvalues: 1D array of eigenvalues
  • eigenvectors: 2D array (column i is eigenvector for eigenvalue i)
  • matrix_type: 'exact' or 'approximate'
  • num_eigenvalues: number computed
  • which_eigenvalues: 'smallest', 'largest', or 'both'
  • Metadata: git_hash, timestamp

• Implementation note: Map user-friendly values to scipy's which parameter:

  • 'smallest' → 'SA' (Smallest Algebraic - most negative eigenvalues)
  • 'largest' → 'LA' (Largest Algebraic - most positive eigenvalues)
  • 'both' → compute smallest and largest separately, concatenate results

• Documentation must emphasize: "smallest" means algebraically smallest (closest to -∞), not smallest magnitude

• Memory management: Store only requested number of eigenvalues/vectors

Testing:

• Unit: Test eigendecomposition on small known matrices (identity, Pauli matrices)
• Unit: Verify orthonormality of eigenvectors
• Unit: Verify eigenvalue ordering
• Unit: Test num_eigenvalues = "all" (full decomposition) on small matrix
• Unit: Test num_eigenvalues = 5 (partial decomposition) on larger matrix
• Unit: Test case insensitivity ("all", "All", "ALL" all work)
• Integration: Test eigendecomposition_matrices='exact' (one output file)
• Integration: Test eigendecomposition_matrices='approximate' (one output file)
• Integration: Test eigendecomposition_matrices='both' (two output files)
• Integration: Test 'smallest', 'largest', and 'both' for which_eigenvalues
• Integration: Verify full vs partial decomposition give same results (for overlap)
• Integration: Verify no unnecessary matrix computations

---

Branch 3: bkk_error_metrics (Error Analysis)

Purpose: Compute three error metrics comparing exact vs approximate

Builds on: bkk_eigendecomposition

Files to modify:

• analysis/analysis.py: Add error_analysis() function
• analysis/config_types.py: Add error analysis fields to AnalysisConfiguration:
  • error_num_eigenvalues: int (how many eigenvalues to compare; 0 = skip eigenvalue errors, default=0)
  • error_matrix_norms: str or list of str (which matrix norms to compute; None = skip matrix errors, default=None)
    • Supported values: "frobenius", "spectral"
    • Can be single string: "frobenius"
    • Can be list: ["frobenius", "spectral"]
  • error_state_inputs: str or list of str (state files for error comparison; None = skip state errors, default=None)
• analysis/README.md: Document error analyses
• analysis/examples/config_full_analysis.py: Add error analysis configuration

New files:

• analysis/tests/test_error_analysis.py: Tests for error metrics

Key implementation details:

Error analysis has three independent components, each explicitly enabled:

1. Eigenvalue Error (enabled when error_num_eigenvalues > 0):

   • Requires both exact and approximate eigendecompositions (will compute if not already done)
   • Compare smallest k eigenvalues from exact vs approximate
   • Output: eigenvalue_errors array with differences exact_eigs[i] - approx_eigs[i]
   • Report both relative and absolute errors
   • Typical use: error_num_eigenvalues = 1 for ground state energy comparison

2. Matrix Norm Error (enabled when error_matrix_norms is not None):

   • Computes specified norm(s) of the difference matrix ||U_exact - U_approximate||
   • Critical scaling challenge: Computing norms requires comparing the two full matrices
   • Two implementation paths:

   Small systems (≤15 qubits): Direct computation

   • Materialize both dense matrices
   • Compute difference: diff = exact - approximate
   • Compute norms directly

   Large systems (>15 qubits): Matrix-free computation (expensive but feasible)

   • Frobenius norm: Can compute without materializing full matrices

     • ||A||_F^2 = sum_ij |A_ij|^2 = trace(A† A)
     • Compute via matrix-free method: apply both operators to each basis state
     • For N qubits: requires 2^N matrix-vector products (expensive but doable)
     • Example: 20 qubits = 1M matvecs, may take minutes to hours
     • Warn user about computation time before starting

   • Spectral norm: Very expensive but computable

     • ||A||_2 = largest singular value = sqrt(largest eigenvalue of A† A)
     • Use power iteration with matrix-free operators
     • Iterative method converges in ~10-100 iterations typically
     • Each iteration requires matrix-vector products with both operators
     • Can take hours for 20+ qubit systems
     • Warn user about computation time before starting

   Implementation approach for Branch 3:

   • Small systems (≤15 qubits): Direct computation with dense matrices (fast)

   • Large systems (>15 qubits): Matrix-free computation with progress warnings

     • Before starting: Log estimated time/cost
     • During computation: Show progress updates
     • Warning message: "Computing Frobenius norm for 20 qubits requires ~1M matrix operations, estimated time: X minutes"

   • User can always Ctrl+C if it's too slow

   • Suggest state-dependent errors as faster alternative in documentation

   • Frobenius norm ("frobenius"): ||A||_F = sqrt(sum of squares of all elements)

   • Spectral norm ("spectral"): ||A||_2 = largest singular value

   • Normalize input to list internally (handle both string and list inputs)

3. State-Dependent Error (enabled when error_state_inputs is not None):

   • Input: string or list of strings specifying initial state files (e.g., "state.npy" or ["state1.npy", "state2.npy"])
   • This is the best-scaling error metric for large systems
   • For each input state:
     • Apply approximate operator: approx_final = approx_operator @ initial_state
     • Apply exact operator: exact_final = exact_operator @ initial_state
     • Compute error: ||exact_final - approx_final||_2
     • Report both absolute and relative to state norm
   • Key scaling advantage:
     • State vectors are O(2^N) not O(2^(2N)) like matrices
     • Can reach 30 qubits (16 GB per state vector)
     • Applies operators without materializing matrices (matrix-free)
   • Implementation:
     • For small systems: use materialized matrices if available (fast)
     • For large systems: use matrix-free operators (PauliStringOperator)
     • Automatically select based on dimension
   • This is independent of numerical_simulation_inputs and exact_simulation_inputs
   • Pattern matches numerical simulation: provide filenames to enable the analysis

Output format (npz or JSON file, e.g., error_analysis.npz):

• eigenvalue_errors: Array of eigenvalue differences (if eigenvalue error enabled)
• eigenvalue_errors_relative: Relative errors (if eigenvalue error enabled)
• matrix_frobenius_error: Scalar (if "frobenius" in error_matrix_norms)
• matrix_spectral_error: Scalar (if "spectral" in error_matrix_norms)
• state_errors: List of dicts with {input_file, absolute_error, relative_error} (if state error enabled)

Memory management:

• Free exact_matrix after computations if no longer needed
• Free unitary_matrix after computations if no longer needed
• Check at end of analyze_algorithm() and explicitly delete large arrays

Testing:

• Unit: Test error calculations with known matrices
• Unit: Test error = 0 when exact = approximate
• Unit: Test error increases with perturbations
• Unit: Test Frobenius norm calculation (both dense and matrix-free paths)
• Unit: Test spectral norm calculation (both dense and matrix-free paths)
• Unit: Test error_matrix_norms as string ("frobenius") and as list (["frobenius", "spectral"])
• Unit: Verify warning messages displayed for large systems
• Integration: Test each error type independently (eigenvalue only, matrix only, state only)
• Integration: Test all three error types together
• Integration: Test matrix norms on system just above dense threshold (16 qubits) - verify matrix-free path works
• Integration: Test error when required data missing (should fail with clear message)
• Integration: Verify errors NOT computed when parameters are None/0
• Performance: Document computation time for matrix norms at 15, 16, 18, 20 qubits

---

Branch 4: bkk_exact_simulation (Complete Pipeline)

Purpose: Numerical simulation using exact matrix

Builds on: bkk_error_metrics

Files to modify:

• analysis/analysis.py: Add exact_numerical_simulation() function
• analysis/config_types.py: Add to AnalysisConfiguration:
  • exact_simulation_inputs: string or list (like numerical_simulation_inputs)
• analysis/README.md: Document exact simulation
• analysis/examples/config_full_analysis.py: Update to demonstrate all analyses together

Key implementation details:

• Very similar to existing numerical_simulation() function
• Uses exact operator (matrix or matrix-free) instead of unitary
• Output files: input_exact_final.npy (distinguish from input_final.npy)
• Can be used independently OR in conjunction with approximate simulation for comparison
• Scales to 25-30 qubits using matrix-free operator
• For large systems:
  • Don't materialize full matrix
  • Apply Hamiltonian via Pauli string evaluation: for each Pauli term, apply to state
  • Memory: state vector only (O(2^N)), not matrix (O(2^(2N)))
  • Same approach used in state-dependent error analysis

Testing:

• Unit: Test exact simulation with identity matrix
• Unit: Test with known Pauli matrix
• Integration: Compare exact simulation with known analytical result
• Integration: Test multiple input states
• End-to-end: Run full analysis with all features enabled

---

Cross-Cutting Concerns

Memory Management Pattern

Key architectural change for scale: Don't always materialize full matrices. Use operators (dense matrix or matrix-free) depending on system size.

def analyze_algorithm(...):
    # Determine system dimension
    num_qubits = hamiltonian.num_qubits()
    dimension = 2 ** num_qubits
    use_dense = (dimension <= 32768)  # 15 qubits threshold
    
    # Track what we need
    needs_approximate = (
        matrix_output_file or 
        numerical_simulation_inputs or 
        error_matrix_norms or
        error_state_inputs or
        (num_eigenvalues > 0 and eigendecomposition_matrices in ['approximate', 'both'])
    )
    needs_exact = (
        exact_matrix_output_file or 
        exact_simulation_inputs or 
        error_matrix_norms or
        error_state_inputs or
        (num_eigenvalues > 0 and eigendecomposition_matrices in ['exact', 'both'])
    )
    
    # Create operators (not always full matrices!)
    if needs_approximate:
        if use_dense:
            approximate_op = _compute_unitary_matrix(...)  # Dense numpy array
        else:
            approximate_op = algorithm  # Algorithm bloq with matvec (matrix-free)
    
    if needs_exact:
        if use_dense:
            exact_op = _compute_exact_matrix_dense(...)  # Dense numpy array
        else:
            exact_op = _create_pauli_operator(hamiltonian)  # PauliStringOperator (matrix-free)
    
    # Use operators for analyses (analyses adapt to dense vs matrix-free)
    if matrix_output_file:
        results["matrix_output"] = output_unitary_matrix(..., approximate_op, use_dense)
    if exact_matrix_output_file:
        results["exact_matrix_output"] = output_exact_matrix(..., exact_op, use_dense)
    if num_eigenvalues > 0:
        results["eigendecomposition"] = eigendecomposition_analysis(..., exact_op, approximate_op)
    if error_state_inputs:
        results["state_errors"] = compute_state_errors(..., exact_op, approximate_op)
    if error_matrix_norms:
        results["matrix_errors"] = compute_matrix_errors(..., exact_op, approximate_op, use_dense)
    
    # Clean up (only needed for materialized matrices)
    if use_dense:
        del approximate_op, exact_op

Design principle: Analyses accept operators that may be dense matrices OR LinearOperators with matvec(). Each analysis adapts its implementation based on operator type.

Configuration Validation

• Each analysis function validates its specific requirements
• analyze_algorithm() validates at least one analysis is requested
• Provide clear error messages for missing dependencies (e.g., error analysis requires both exact and unitary)

CI/CD Updates

Both .github/workflows/python-app.yml and .gitlab-ci.yml:

• Already have run-analysis-tests jobs that run pytest
• New tests automatically included via pytest tests/
• No structural changes needed
• Verify each branch passes before merge

Documentation Strategy

README.md sections to add/update:

• "Matrix Output" section: Add exact matrix subsection
• New "Eigendecomposition Analysis" section
  • Explain num_eigenvalues options:
    • Small positive int (e.g., 5): for large systems, compute only a few eigenvalues (RECOMMENDED for production)
    • "all": compute full eigendecomposition (all eigenvalues/eigenvectors) for small systems only
    • 0: disabled
  • Example usage: analysis.num_eigenvalues = "all" or analysis.num_eigenvalues = 5
  • Must clearly explain: 'smallest' = algebraically smallest (most negative, ground state energy)
  • Example: For Hamiltonian with eigenvalues [-10, -5, 0, 5, 10], 'smallest' gives [-10, -5, ...]
  • Clarify this is NOT smallest magnitude (which would give [0, 5, -5, ...])
  • System size guidance:
    • ≤10 qubits: "all" is feasible
    • 10-12 qubits: "all" possible but expensive
    • ≥12 qubits: use partial (specify small k)
    • 20 qubits (production scale): MUST use partial, "all" not feasible
    • Partial eigendecomposition scales well beyond 20 qubits
• New "Error Analysis" section with subsections for each error type
  • Clearly document "frobenius" vs "spectral" norms:
    • Frobenius: element-wise difference, fast to compute
    • Spectral: worst-case effect on any quantum state, more physically meaningful but expensive
    • Include when to use each (Frobenius for quick comparisons, spectral for physical bounds)
• New "Exact Numerical Simulation" section
• Update "Configuration Options" summary

Code comments:

• Docstrings for all new functions
• Inline comments for non-obvious algorithms (e.g., sparse matrix construction)
• NO intermediate process documentation in repo

Example updates:

• config_full_analysis.py (renamed in Branch 0) evolves across branches to demonstrate cumulative features
• Each branch adds its new analysis capabilities to this example
• Final version demonstrates complete analysis pipeline with all available features

Testing Strategy

Unit tests:

• Each new module gets comprehensive unit tests
• Test edge cases: empty matrices, single-qubit, identity, etc.
• Test error handling: dimension mismatches, missing data

Integration tests:

• Each branch adds integration test for its new capability
• Integration tests use small systems (2-4 qubits) for CI/CD speed
• At least one test exercises full pipeline end-to-end

Production scale expectations:

• Normal production use: 20 qubits (dimension 2^20 = 1,048,576)
• Target: 25-30 qubits (dimension up to 2^30 = 1 billion)
• Stretch goal: Push toward 50 qubits where feasible
• Higher is always better - users will push to limits

Scaling limits by operation type:

• Dense matrix storage: Limited to ~15 qubits (16 GB at 15 qubits, 16 TB at 20 qubits)
• State vectors: Can reach ~30 qubits (16 GB for 30-qubit state)
• Sparse/matrix-free operations: Can reach 25-30+ qubits depending on Hamiltonian sparsity

Design principles for maximum scale:

• Avoid materializing full dense matrices for large systems
• Matrix-free operators: Compute matrix-vector products directly from Pauli strings
• Sparse eigensolvers: Use iterative methods that only need matrix-vector products
• User control: Let user choose dense vs sparse representations based on system size
• Graceful degradation: Operations fail cleanly with clear messages when limits exceeded

Test organization:

• tests/test_matrix_operations.py: Matrix utilities
• tests/test_exact_matrix.py: Exact matrix computation
• tests/test_eigendecomposition.py: Eigenanalysis
• tests/test_error_analysis.py: Error metrics
• tests/test_analysis_integration.py: Update to include new features in full pipeline test

Dependencies Between Branches

main (bkk_numerical)
    ↓
bkk_config_full_analysis (rename example, establish baseline)
    ↓
bkk_exact_matrix (compute exact H matrix)
    ↓
bkk_eigendecomposition (eigenanalysis of both matrices)
    ↓
bkk_error_metrics (compare exact vs approximate)
    ↓
bkk_exact_simulation (simulate with exact matrix)

Each branch:

1. Builds on previous work
2. Can be merged independently as a complete feature
3. Has its own tests and documentation
4. Results in a working, testable system

Merge Request Checklist (for each branch)

• All unit tests pass locally
• All integration tests pass locally
• GitHub Actions CI passes
• GitLab CI passes
• README.md updated for new features
• Example config updated
• Code review requested
• Dependencies on previous branches clearly stated in MR description
• Memory management verified (no leaks in large systems)

Progress Tracking File

Location: /Users/bkkrueger/.claude/plans/exact_matrix_implementation_progress.md

This file will track:

• Current branch being worked on
• Completed branches and their merge commit hashes
• Any issues encountered and resolutions
• Notes for future work
• Branch dependencies and merge order

Format:

# Exact Matrix Implementation Progress

## Status: [In Progress / Paused / Complete]
Current branch: bkk_exact_matrix
Last updated: YYYY-MM-DD

## Completed Branches
- None yet

## In Progress: bkk_exact_matrix
- [x] Created branch from bkk_numerical
- [x] Added matrix_operations.py
- [ ] Unit tests for matrix operations
- [ ] Integration tests
- [ ] Documentation updates
- [ ] CI/CD verification

## Next Steps
1. Complete bkk_exact_matrix
2. Create MR for bkk_exact_matrix
3. Begin bkk_eigendecomposition

## Notes
- Remember to test with both small (2-3 qubit) and medium (4-6 qubit) systems
- Watch memory usage for large matrices
- Coordinate with team on merge timing

Critical Files Reference

Starting point: /Users/bkkrueger/research/quantum_computing/lisdi-qre/qhat/worktrees/bkk_numerical/

Key files to understand:

• analysis/analysis.py:237-291: analyze_algorithm() main dispatcher
• analysis/analysis.py:64-94: _compute_unitary_matrix() pattern to follow
• analysis/analysis.py:152-234: numerical_simulation() pattern for exact simulation
• analysis/config_types.py:157-168: AnalysisConfiguration class
• analysis/file_io.py: save_matrix, load_state patterns to reuse
• analysis/hamiltonian.py:132-169: get_all_pauli_strings() method
• hamiltonian_generator/eigendecompose.py: Existing eigendecomposition code (for reference)
• common/logging_utils.py: Centralized logging setup (new in main)
• analysis/tests/conftest.py: Test logging configuration (new in main)

Logging patterns (after merging main into bkk_numerical):

# At top of module
import logging
logger = logging.getLogger(__name__)

# In functions
logger.info("Starting computation...")  # replaces config_general.log()
logger.verbose("Detailed info...")      # replaces config_general.log_verbose()
logger.debug("Debug info...")           # replaces config_general.log_debug()

CI/CD:

• .github/workflows/python-app.yml:80-90: analysis unit tests job
• .gitlab-ci.yml:83-94: analysis tests job

Estimated Complexity

Branch 0 (config baseline): Very Low

• Simple rename and documentation update
• No new functionality, just organization
• Should be quick to complete

Branch 1 (exact matrix): Medium

• Core matrix operations are straightforward
• Testing with sparse matrices adds complexity
• Foundation for other branches

Branch 2 (eigendecomposition): Low-Medium

• scipy provides robust eigenvalue solvers
• Main complexity is handling different request modes

Branch 3 (error metrics): Medium

• Three different error types to implement
• State-dependent error requires careful coordination
• Memory management most critical here

Branch 4 (exact simulation): Low

• Largely copy of existing numerical_simulation()
• Integration and testing is main work

Scaling Summary by Analysis Type

Analysis	Small (≤15 qubits)	Medium (16-20 qubits)	Large (21-30 qubits)	Notes
Exact matrix output	✓ Dense	✗ Too large	✗ Too large	Limited by dense matrix storage (16 GB at 15 qubits)
Approx matrix output	✓ Dense	✗ Too large	✗ Too large	Same limit as exact matrix
Eigendecomposition (full)	✓ Dense	✗ Too large	✗ Too large	Requires dense matrix
Eigendecomposition (partial)	✓ Dense or sparse	✓ Sparse	✓ Sparse	Matrix-free operator, scales well
Exact simulation	✓ Dense	✓ Matrix-free	✓ Matrix-free	State vector limit ~30 qubits (16 GB)
Approx simulation	✓ Dense	✓ Algorithm op	✓ Algorithm op	Same as exact simulation
Eigenvalue error	✓ Dense	✓ Sparse	✓ Sparse	Uses partial eigendecomposition
Matrix norm error	✓ Dense	⚠ Slow (matrix-free)	⚠ Very slow (matrix-free)	Frobenius: 2^N ops; Spectral: iterative
State error	✓ Dense	✓ Matrix-free	✓ Matrix-free	Best-scaling error metric

Key takeaways:

• Dense operations limited to ~15 qubits (matrix output, full eigendecomp)
• Matrix norm errors work but slow for >15 qubits (minutes to hours for 20+ qubits)
• Partial eigendecomposition reaches 25-30 qubits via matrix-free operators
• Simulations and state errors reach 25-30 qubits (limited by state vector size)
• Production recommendation for 20+ qubits:
  • Use partial eigendecomposition for eigenvalue comparison
  • Use state errors for best performance
  • Use matrix norms if needed, but expect long computation times

Success Criteria

Each branch:

• Passes all existing tests (no regressions)
• Passes new tests (>90% coverage for new code)
• Documentation complete
• CI/CD green on both platforms

Overall project:

• User can run full analysis comparing exact vs approximate
• Example config demonstrates all capabilities
• All error metrics validated on test systems
• No memory leaks with realistic problem sizes (20+ qubits)
• Performance validated at production scale:
  • Confirm partial eigendecomposition works smoothly at 20 qubits
  • Verify memory usage is reasonable for 20-qubit systems
  • Document any system-size limitations discovered
• Ready for scientific use and publication