README: Tensor Network HHL Simulation

Project Title

Solving Systems of Linear Equations: HHL from a Tensor Networks Perspective

This repository accompanies the work presented in the paper:

"Solving Systems of Linear Equations: HHL from a Tensor Networks Perspective" Alejandro Mata Ali and Iñigo Perez Delgado and Marina Ristol Roura and Aitor Moreno Fdez. de Leceta and Sebastián V. Romero (2025) arXiv:2309.05290

It can be consulted online in the Streamlit webpage: https://tensornetworks-hhl-algorithm.streamlit.app/

The project implements a classical simulation of the quantum Harrow-Hassidim-Lloyd (HHL) algorithm using tensor networks and qudit formalism. The goal is to provide a quantum-inspired solver that models the ideal behavior of HHL efficiently on classical hardware, enabling benchmarking and theoretical lower-bound estimations.

---

Files

• tensor_network_HHL.ipynb: Jupyter notebook containing all code, explanations, experiments, and plots. It reproduces the results presented in the paper.

---

Requirements

Install the dependencies with:

pip install numpy matplotlib scipy torch qiskit qiskit_ibm_runtime qiskit_aer

The notebook is compatible with standard Python 3.x and requires no GPU. All computations were tested on CPU.

---

Usage

Open the notebook in Jupyter:

jupyter notebook tensor_network_HHL.ipynb

You may execute all cells sequentially to:

1. Define the tensor operations for the TN-HHL algorithm.
2. Construct tensors for QPE, inversion, and evolution operators.
3. Apply the method to benchmark problems:

• Forced harmonic oscillator
• Forced damped oscillator
• 2D static heat equation with sources

4. Compare TN-HHL performance to:

• Exact inversion (PyTorch)
• Qiskit HHL simulation (for small cases)

Each section is self-contained and annotated for clarity.

---

Summary of the Algorithm

• The notebook encodes the HHL quantum circuit using tensor networks.
• It implements all gates (QPE, inversion, unitaries) as tensor contractions.
• The final solution vector $\vec{x}$ is obtained deterministically, bypassing quantum limitations like post-selection.

Parameter conventions

• tau (API name t) sets the spectral grid spacing: $\Delta\lambda = 1/\tau$.
• mu (API names n_eigen / num_eigen) is the phase-register dimension. It is not by itself the spectral resolution.
• The signed non-aliased spectral range depends on both parameters: $|\lambda| < \mu/(2\tau)$, equivalently $|\tau\lambda| < \mu/2$.
• n_c is used only for a binary phase register with $\mu = 2^{n_c}$. The general tensor implementation does not require $\mu$ to be a power of two.
• In the quantum-circuit HHL, the adaptive experiment rule is $C_{\mathrm{phys}}(\tau)=\eta\min(\lambda_{\min},1/\tau)$ with $0<\eta\leq1$, so $C_{\mathrm{bin}}=\tau C_{\mathrm{phys}}\leq\eta$. The general rotation function still rejects arbitrary constants outside the arcsin domain without clipping.

Fourier and inversion conventions

• The paper derivation uses the unnormalized Fourier matrix $H[a,b]=\exp(2\pi iab/\mu)$, with $H^{-1}=H^\dagger/\mu$.
• The code uses the normalized QFT $F=H/\sqrt{\mu}$, which is unitary.
• The inverter tensor stores the global factor $\tau/\mu$ so that the normalized implementation remains equivalent to the paper's unnormalized formulation.

Real-valued experimental scope

The numerical experiments accompanying the current paper use real-valued matrices and right-hand sides. The present public implementation returns real-valued solution vectors and has been validated only for this real-valued benchmark setting. The underlying tensor-network formulation is not restricted to real systems. Supporting genuinely complex-valued inputs would require retaining the complete complex output and using phase-aware validation metrics.

---

Reproducibility

The notebook contains the explanatory examples and the original interactive workflow. The additional experiments requested by Reviewer 1 are generated by a separate, reproducible Python command and then loaded by the final notebook cells.

Reproducing the reviewer-1 spectral and discretization diagnostics

Run the following commands from the project root. On the first execution in Windows PowerShell, activate the virtual environment and install the additional dependencies:

.\.venv\Scripts\Activate.ps1
python -m pip install -r requirements.txt -r requirements_experiments.txt

On Linux or macOS, use:

source .venv/bin/activate
python -m pip install -r requirements.txt
python -m pip install -r requirements_experiments.txt

After installation, generate all Reviewer 1 experiments with:

python -m experiments.run_reviewer_r1_c5_c6

This command normally takes approximately 30 seconds to 2 minutes, depending mainly on available memory and CPU performance. It recreates the scientific tables and six PDF figures in artifacts/reviewer_r1_c5_c6/: application results, spectral diagnostics, all values of tau*lambda_j, parameter-search data, random-matrix results, the sensitivity sweep, and TN/filter validation.

Wait for the message Generated scientific results ... before opening tensor_network_HHL.ipynb and selecting Run All. The order is therefore:

1. Run python -m experiments.run_reviewer_r1_c5_c6.
2. Wait until result generation finishes.
3. Open the notebook and run all cells.

The notebook's final sensitivity-analysis cells read artifacts/reviewer_r1_c5_c6/hyperparameter_sweep.csv, so running the command first is required. Running the notebook is not required to generate the Reviewer 1 results; it provides the interactive explanations and visualizations. Conversely, the reviewer command does not execute the notebook's original examples or its controlled Qiskit comparison.

The corrected 20-instance Qiskit-versus-TN comparison is generated independently of the notebook with

python -m experiments.run_reviewer_r1_c7

Both implementations use the same selected (mu, tau) and the same finite evolution

U = exp(2*pi*i*tau*A/mu).

The exact statevector produces two distinct postselections. Conditioning only on ancilla success and tracing out the phase register gives the generally mixed density matrix $\rho_{\mathrm{sys}\mid a=1}$ with probability $p_{a=1}$. Projecting jointly onto ancilla success and clock zero gives the pure state $|x_{a=1,c=0}\rangle$ with probability $p_{a=1,c=0}$. The primary TN comparison is the squared overlap with this pure joint branch. Fidelity against $\rho_{\mathrm{sys}\mid a=1}$, its purity, and probability RMSE are secondary diagnostics. The separate seeded 100,000-shot run is conditioned only on the ancilla and is not presented as a direct TN-state verification.

For this comparison the selector uses mu = 128, 256, 512, 1024, 2048, 40 logarithmic tau values per mu, and $\eta=0.9$. It rejects aliasing and insufficient separation from the singular bin and validates every rotation before circuit construction. All 20 selected configurations have right-hand-side relative filter error at most 1%. For fixed (mu, tau), changing $C$ does not change the normalized filter state; it changes only rotation feasibility and success probability. The TN-equivalent joint probability includes StatePreparation's normalization: $C_{\mathrm{phys}}^2\lVert f_{\mu,\tau}(A)b\rVert^2/\lVert b\rVert^2$.

Reviewer 2 comments 3 and 4 add isolated-process memory measurements and a limited empirical scaling study:

python -m experiments.run_reviewer_r2_c3_c4 --refresh-comparison-memory

The flag explicitly remeasures isolated-process RSS for all 20 Qiskit–TN rows before running the scaling study. Without it, the command requires the comparison CSV to contain complete memory columns. Baseline RSS, absolute peak RSS, and their difference are retained, but the primary empirical memory metric in summaries and figures is the peak increment $\Delta\mathrm{RSS}=\mathrm{RSS}{\mathrm{peak}}-\mathrm{RSS}{\mathrm{baseline}}$.

The default $N$ sweep uses N = 16, 32, 64, 96, 128, 192, 256 at mu=64. The default phase-dimension sweep uses mu = 16, 32, 64, 128, 256, 512, 1024, 2048 at N=32; both use tau=20, five timing repetitions, and five isolated memory repetitions. Alternative values can be passed with --n-values, --mu-values, --n-sweep-mu, and --mu-sweep-n.

The public implementation explicitly materializes inverse_phase_kickback with shape $\mu\times N\times N$. Its dense preparation performs repeated matrix–matrix products and therefore includes $O(\mu N^3)$ work. The runner also reports dominant_tensor_storage_estimate_bytes, an algebraic estimate for explicitly materialized scientific tensors with dominant storage $O(\mu N^2+\mu^2+N^2)$. This estimate is not measured RSS and excludes temporary allocations internal to LAPACK, BLAS, PyTorch, and matrix_exp. The fitted log-log exponents are descriptive observations over the finite measured range, not proofs of asymptotic order. The independently timed spectral filter remains a validation reference.

The oscillator discretizations use 100 intervals and 101 nodes; the two prescribed endpoints leave 99 interior unknowns. The damped system therefore has a 99 by 99 physical matrix and a 198 by 198 Hermitian embedding.

The (mu, tau) pairs are selected from a versioned, predefined logarithmic grid, now extending through mu=4096, using the predicted right-hand-side error of the HHL spectral filter. The selector is labeled exact_spectral_benchmark_oracle: it uses a full eigendecomposition and the classical reference solution as an offline calibration tool for controlled experiments. It is not presented as a practical parameter selector for an unknown large system. Its measured time is recorded separately from TN and reference-solver timings.

Random matrices use the fixed seed 12345. Sensitivity figures show two-sided 95% Student-t confidence intervals for the mean across 20 instances. The full 81-by-20 sensitivity sweep uses the algebraically equivalent spectral filter and is labeled spectral_filter_equivalent; it is not described as 1,620 TN contractions. A configured representative subset is executed with the complete TN contraction and written to tn_filter_validation.csv, which records TN/filter discrepancies and separate timings. Selected application and random runs also use the complete TN contraction.

The benchmark-oracle selection cost is approximately O(N^3 + sum_c(N*mu_c + N^2)) for the evaluated candidates, including the eigendecomposition. This calibration cost is distinct from the fixed-parameter TN contraction cost and must be included separately in any end-to-end runtime discussion. See [scientific_addings.md](scientific_addings.md) for the scientific interpretation and manuscript-facing caveats.

The versioned Python modules and the reviewer command above are the authoritative source for these additional diagnostics.

---

Reference

If you use this code, please cite the original paper (also encoded in [CITATION.cff](CITATION.cff)):

@misc{ali2024solvingsystemslinearequations,
      title={Solving Systems of Linear Equations: HHL from a Tensor Networks Perspective}, 
      author={Alejandro Mata Ali and Iñigo Perez Delgado and Marina Ristol Roura and Aitor Moreno Fdez. de Leceta and Sebastián V. Romero},
      year={2024},
      eprint={2309.05290},
      archivePrefix={arXiv},
      primaryClass={quant-ph},
      url={https://arxiv.org/abs/2309.05290}, 
}

---

License and copyright

This software is released under the MIT License. Copyright (c) 2025 Alejandro Mata Ali — see [LICENSE](LICENSE).

The accompanying paper (arXiv:2309.05290) lists additional coauthors (Iñigo Perez Delgado, Marina Ristol Roura, Aitor Moreno Fdez. de Leceta, and Sebastián V. Romero). Paper authorship is for scientific credit and citation; it does not change the software copyright holder named in LICENSE.

Third-party runtime dependencies (including Qiskit and PyTorch) are described in [NOTICE](NOTICE).