PySpec

PySpec reconstructs smeared spectral densities from lattice correlators. Current version supports the Hansen-Lupo-Tantalo (HLT) method. Given a target smearing kernel K(omega) it computes the coefficients g_t which provides the best approximation of that kernel by balancing systematic and statistical error.

Everything is evaluated in arbitrary precision with mpmath.

Distribution: PySpec  ·  importable as pyspec.

---

Table of Contents

• The method
• Features
• Requirements
• Installation
• Quickstart
• Kernels
• Constraints
• Tuning lambda
• Project layout
• License
• Bug reports
• Author and contact

---

The method

The lattice correlator is a Laplace transform of the spectral density,

C(t) = int_{E0}^{inf} b_t(omega) rho(omega) d omega,
   b_t(omega) = exp(-(t+1) omega)   [ + exp(-(T-t-1) omega) with periodic bc ]

Following arXiv:1903.06476 PySpec looks for coefficients g such that sum_t g_t b_t(omega) reproduces a chosen smearing kernel K(omega), so that sum_t g_t C(t) is an estimator of the smeared density. The coefficients minimize

(1 - lambda) A[g]/A[0] + lambda B[g]

where A[g] measures the distance between the reconstructed and the target kernel and B[g] = g^T Cov g is the variance of the result. The minimum has a closed form in the presence of linear constraints V g = c.

Features

Reconstruction

• Reconstructor: one instance for each kernel and set of lattice parameters. The kernel-independent block A_{tr} is computed once, and the LU factorization of W = (1 - lambda) A + lambda Cov is cached per lambda, so scanning lambda or swapping kernel is cheap.
• Periodic and non-periodic boundary conditions.
• set_kernel recomputes only the kernel-dependent terms, for scans over the energy at which the density is evaluated.

Kernels

• Currently supported kernels: SigmoidKernel, MomentsKernel, ExpKernel, GaussDKernel, ExpMomentsKernel, SechKernel, ExpSigmoidKernel, ExpSigmoidModKernel, AsymSechKernel.
• Custom kernels by subclassing Kernel and implementing evaluate: the Laplace transform, the square integral and the area are then obtained by quadrature, and can be overridden when a closed form is known.

Constraints

• NormalizationConstraint: the reconstructed kernel has the same area as the target one.
• PointMatchConstraint: the two kernels agree exactly at a chosen energy.
• LinearConstraint: template for any constraint of the form v^T g = c.

Tuning

• Automatic selection of lambda by comparing reconstructions obtained with different constraint sets, with a plateau search and a compatibility threshold in units of the statistical error.

---

Requirements

Python >= 3.13, with mpmath, numpy and scipy: installed automatically.

---

Installation

pip install "git+https://github.com/laudid46/pyspec.git@v1.0"

---

1. Clone the repository

git clone https://github.com/laudid46/pyspec.git
cd pyspec

---

2. Create and activate a virtual environment (recommended)

python3.13 -m venv .venv
source .venv/bin/activate           # macOS / Linux
# .\.venv\Scripts\Activate.ps1      # Windows PowerShell

---

3. Install the package

For regular/development (changes to the source are picked up immediately):

pip install .                # regular installation
pip install -e .             # editable installation
pip install -e ".[test]"     # editable, with the test dependencies

---

4. Check the installation

From the repository root:

pytest               # test package installation

This requires the test extra (see above). Expected output: 37 passed.

---

Quickstart

import numpy as np
import mpmath as mp
from pyspec import Reconstructor, GaussDKernel, NormalizationConstraint, PointMatchConstraint

mp.mp.dps = 60                       # the inverse problem needs the digits

C   = ...                            # correlator, C[t] for t = 0, 1, ...
Cov = ...                            # its covariance matrix, NOT normalized

kernel = GaussDKernel(omega_star=1.0, sigma=0.3)

rec = Reconstructor(kernel, [C, Cov],
                    tmin=0, tmax=10,      # timeslices used
                    E0=0.05,              # lower bound of integrations
                    T=None,               # lattice time extent, None if not periodic
                    cov_norm=float(C[0]**2))

# coefficients at a fixed lambda
g = rec.coefficients(mp.mpf('0.5'), [NormalizationConstraint()])

# the reconstructed kernel and the smeared density
rec.K(1.0, mp.mpf('0.5'), [NormalizationConstraint()])
rec.rho(mp.mpf('0.5'), [NormalizationConstraint()])

# or let PySpec tune lambda for you
out = rec.run([[NormalizationConstraint()],
               [NormalizationConstraint(), PointMatchConstraint(1.2)]],
              tuning=True, N=32)

print(out['opt_info']['rho'], '+/-', out['opt_info']['sigma'])
print('lambda =', out['lambda_opt'])

tmin and tmax select the timeslices t = tmin, ..., tmax-1, which use the correlator entries C[tmin+1], ..., C[tmax]: the t = 0 entry is the normalization, not a basis function.

---

Kernels

Every kernel exposes four primitives in arbitrary precision:

method	meaning
evaluate(omega)	K(omega)
laplace(s, E0)	int_{E0}^{inf} exp(-s omega) K(omega) d omega
square_integral(E0)	int_{E0}^{inf} K(omega)^2 d omega
integral(E0)	int_{E0}^{inf} K(omega) d omega, i.e. laplace(0, E0)

A new kernel only needs evaluate; the rest falls back to adaptive quadrature. If the kernel has a cusp, override laplace and square_integral and split the integration path on it, otherwise the quadrature silently loses several digits.

---

Constraints

A linear constraint is a pair (v, c) enforcing v^T g = c. To add one, subclass LinearConstraint and implement get_v, get_value and the key property, which must return a string identifying the constraint uniquely: it is used as the cache key.

from pyspec import LinearConstraint

class MyConstraint(LinearConstraint):
    def get_v(self, kernel, tmin, tmax, E0, T=None):
        ...                                  # mp.matrix of shape (tmax-tmin, 1)
    def get_value(self, kernel, E0):
        ...                                  # mp.mpf
    @property
    def key(self):
        return 'MyConstraint'

---

Tuning lambda

run scans lambda, converts each point into the corresponding A[g]/A[0], looks for the plateau of rho, and picks the largest A[g]/A[0] at which two different reconstruction strategies still agree within threshold sigmas. The reconstruction is then re-run at the selected lambda, and the result is compared with the interpolated one.

The returned dictionary contains Aopt, lambda_opt, rho_interp, sigma_interp, the interpolation grid, the full scan for each constraint set, the index of the selected set, and opt_info with the final rho and sigma. See the docstring of run for every option.

---

Project layout

PySpec/
├── LICENSE                    # GPL-3.0
├── README.md                  # this file
├── pyproject.toml             # metadata and dependencies
├── CHANGELOG.md               # version history
├── MANIFEST.in                # files shipped in the source distribution
├── src/pyspec
│       ├── __init__.py        # re-exports the whole public API
│       ├── __version__.py
│       ├── kernels.py         # target smearing kernels K(omega)
│       ├── functionals.py     # basis b_t, R_t, f_t, A, W and the A, B functionals 
│       ├── solver.py          # linear constraints and the closed-form solver 
│       ├── reconstruction.py  # the Reconstructor class and the tuning of lambda
│       └── utils.py           # conversion to mpmath matrices
├── tests/
│   └── test_00.py             # numerical test suite
└── doc/
    └── pyspec_logo.png        # pyspec logo

---

License

Distributed under the GNU General Public License v3.0 (GPL-3.0). See the LICENSE file for the full terms.

---

Bug reports

https://github.com/laudid46/pyspec/issues

---

Author and contact

Davide Laudicina — davide.laudicina1@gmail.com

Repository: https://github.com/laudid46/pyspec