SpectrumInversion¶

Overview¶

The rheojax.transforms.SpectrumInversion transform recovers the continuous relaxation spectrum \(H(\tau)\) from measured dynamic moduli \(G'(\omega)\), \(G''(\omega)\) or relaxation modulus \(G(t)\). This is a classical ill-posed inverse problem in rheology, requiring regularization for stable solutions.

Key Capabilities:

Tikhonov regularization: Automatic parameter selection via GCV (Generalized Cross-Validation) or manual control
Maximum entropy method: Information-theoretic approach preserving positivity
Dual-source input: Works from oscillation data \(G^*(\omega)\) or relaxation data \(G(t)\)
Non-negative spectrum: Physical constraint \(H(\tau) \ge 0\) enforced

The continuous relaxation spectrum provides a material fingerprint that reveals the distribution of relaxation mechanisms—from single-mode Maxwellian fluids (delta function) to broad distributions in polymer melts and filled systems.

Mathematical Theory¶

Continuous Spectrum Definition¶

The continuous relaxation spectrum \(H(\tau)\) is defined through the integral representation:

\[G(t) = G_e + \int_{-\infty}^{\infty} H(\tau) \exp(-t/\tau) \, d(\ln \tau)\]

The dynamic moduli are related to \(H(\tau)\) by the kernel integrals:

\[G'(\omega) = G_e + \int_{-\infty}^{\infty} H(\tau) \frac{\omega^2 \tau^2}{1 + \omega^2 \tau^2} \, d(\ln \tau)\]

\[G''(\omega) = \int_{-\infty}^{\infty} H(\tau) \frac{\omega \tau}{1 + \omega^2 \tau^2} \, d(\ln \tau)\]

Physical interpretation: \(H(\tau) d(\ln \tau)\) is the modulus contribution from relaxation mechanisms with times between \(\tau\) and \(\tau + d\tau\).

The Inverse Problem¶

Recovering \(H(\tau)\) from data is an ill-posed Fredholm integral equation of the first kind. Discretizing on a log-spaced \(\tau\) grid gives:

\[\mathbf{b} = A \mathbf{H} + \boldsymbol{\varepsilon}\]

where \(A\) is the kernel matrix, \(\mathbf{H} = [H(\tau_1), \ldots, H(\tau_M)]\), and \(\boldsymbol{\varepsilon}\) is noise.

Why ill-posed? The singular values of \(A\) decay rapidly—small noise in \(\mathbf{b}\) maps to large oscillations in \(\mathbf{H}\) unless regularized.

Tikhonov Regularization¶

Tikhonov regularization stabilizes the inversion by penalizing solution norm:

\[\min_{\mathbf{H}} \left\{ \|A \mathbf{H} - \mathbf{b}\|^2 + \lambda^2 \|L \mathbf{H}\|^2 \right\}\]

\(\lambda\) controls the trade-off between data fidelity and smoothness
\(L = I\) (zeroth order): penalizes large \(H\) values
\(L = D_1\) (first order): penalizes gradient (promotes smoothness)

Solution: \(\mathbf{H} = (A^T A + \lambda^2 L^T L)^{-1} A^T \mathbf{b}\)

GCV Selection of λ (Implementation Detail)¶

The Generalized Cross-Validation criterion selects \(\lambda\) by minimizing the predicted leave-one-out error without actually performing \(n\) separate inversions:

\[\text{GCV}(\lambda) = \frac{\|A \mathbf{H}_\lambda - \mathbf{b}\|^2} {\left[\text{tr}(I - M(\lambda))\right]^2}\]

where \(M(\lambda) = A(A^T A + \lambda^2 L^T L)^{-1} A^T\) is the influence matrix (also called the hat matrix). The numerator is the residual sum of squares; the denominator penalizes overfitting by measuring the effective number of free parameters.

Physical intuition: GCV balances two extremes:

\(\lambda \to 0\): Zero regularization. \(\mathbf{H}\) fits the noise. Residual is small, but the denominator collapses (all data points are “fit”), so GCV → ∞.
\(\lambda \to \infty\): Infinite regularization. \(\mathbf{H} \to 0\). Residual is large. GCV → ∞ again.
Optimal λ: The minimum of the GCV curve sits between these extremes—enough regularization to suppress noise amplification, but not so much that physical features are lost.

Fast SVD path (L = I): When the regularization matrix is the identity (zeroth-order Tikhonov, which is the RheoJAX default), the implementation uses the SVD of \(A\):

\[A = U \Sigma V^T\]

to express all quantities analytically via the filter factors \(f_j = \sigma_j^2 / (\sigma_j^2 + \lambda^2)\):

\[\|\text{residual}\|^2 = \sum_j (1 - f_j)^2 (U^T \mathbf{b})_j^2 + \|\mathbf{b}_\perp\|^2\]

\[\text{tr}(I - M) = n - \sum_j f_j\]

where \(\mathbf{b}_\perp = (I - U U^T)\mathbf{b}\) is the component of the data orthogonal to the range of \(A\) (constant across all \(\lambda\)). The SVD is computed once at \(O(\min(n, m)^2 \max(n, m))\) cost, after which each of the 50 candidate \(\lambda\) values (log-spaced from \(10^{-6}\) to \(10^4\)) is evaluated in \(O(\min(n, m))\) time—a significant speedup over the general case.

General case (L ≠ I): When a non-identity regularization matrix is used (e.g., first-derivative operator for smoothness), the fast SVD path is not available. Instead, the implementation solves the normal equations for each candidate \(\lambda\) via:

\[\mathbf{H}_\lambda = (A^T A + \lambda^2 L^T L)^{-1} A^T \mathbf{b}\]

and computes the influence trace as \(\text{tr}((A^T A + \lambda^2 L^T L)^{-1} A^T A)\). This costs \(O(m^3)\) per candidate but avoids forming the full \(n \times n\) influence matrix.

Diagnostic: The returned regularization_param in SpectrumResult allows you to verify the selected \(\lambda\). If the spectrum looks over-smoothed, you can override with a smaller manual value; if it oscillates, use a larger one.

Maximum Entropy Method¶

The maximum entropy approach maximizes the Shannon entropy of \(H(\tau)\) subject to a data fidelity constraint:

\[\max_{\mathbf{H}} \left\{ S = -\sum_i H_i \ln(H_i / m_i) \right\} \quad \text{subject to} \quad \chi^2 \le \chi^2_{\text{target}}\]

where \(m_i\) is a prior model (default: uniform). The solution is found via iterative multiplicative updates (Bryan’s algorithm):

\[H_i^{(k+1)} = H_i^{(k)} \exp\!\left( -\lambda \frac{\partial \chi^2}{\partial H_i} \Big/ \left(1 + \lambda \left|\frac{\partial \chi^2}{\partial H_i}\right|\right) \right)\]

Advantages over Tikhonov: Automatically enforces \(H(\tau) > 0\), produces maximally non-committal spectra, and avoids oscillatory artifacts.

Parameters¶

SpectrumInversion Parameters¶
Parameter	Type	Default	Description
`method`	str	`"tikhonov"`	Inversion method: `"tikhonov"` or `"max_entropy"`.
`n_tau`	int	`100`	Number of \(\tau\) points in the output spectrum.
`tau_range`	tuple \| None	`None`	Explicit \((\tau_{\min}, \tau_{\max})\). Auto-detected if `None`.
`regularization`	float \| None	`None`	Manual \(\lambda\). If `None`, auto-selected via GCV.
`source`	str	`"oscillation"`	Input data type: `"oscillation"` or `"relaxation"`.
`G_e`	float	`0.0`	Equilibrium modulus (Pa). Set to 0 for liquids.

Parameter Selection Guidelines¶

Method choice:

Tikhonov: Default, fast, well-understood. Best for clean data with moderate noise.
Maximum entropy: Better for noisy data and when positivity is critical. Slightly slower due to iterative solution.

n_tau:

50–100: Standard (adequate for most polymer systems)
200+: High resolution (for resolving closely-spaced relaxation modes)
<50: Coarse (fast screening)

tau_range auto-detection logic:

When tau_range=None (default), the \(\tau\) grid boundaries are inferred from the input data with a 1-decade safety margin on each side:

Oscillation source (source="oscillation"):

The relationship \(\omega \sim 1/\tau\) maps the frequency window to the relaxation time window with an inversion:

\[\tau_{\min} = \frac{1}{10 \, \omega_{\max}}, \qquad \tau_{\max} = \frac{10}{\omega_{\min}}\]

For example, a frequency sweep from \(\omega = 0.01\) to \(100\) rad/s produces \(\tau \in [10^{-3}, 10^{3}]\) s—six decades, extending one decade beyond the data on each side.

Relaxation source (source="relaxation"):

The time axis maps directly to relaxation times:

\[\tau_{\min} = \frac{t_{\min}}{10}, \qquad \tau_{\max} = 10 \, t_{\max}\]

For example, relaxation data from \(t = 0.001\) to \(100\) s produces \(\tau \in [10^{-4}, 10^{3}]\) s.

Why the 10× safety margin? The kernel functions (\(\omega^2\tau^2/(1+\omega^2\tau^2)\) for \(G'\), \(\omega\tau/(1+\omega^2\tau^2)\) for \(G''\)) have significant sensitivity beyond the strict \(1/\omega\) boundary. A mode at \(\tau = 0.1/\omega_{\max}\) still contributes ~1% to \(G'(\omega_{\max})\), which is resolvable with good SNR. Without the margin, edge modes are truncated and the spectrum appears artificially narrow.

When to override: Use explicit tau_range when:

Data is noisy at the frequency extremes (tighten range to avoid fitting noise)
You know the material has no modes outside a specific range (e.g., single-peak gel)
You want to zoom into a specific part of the spectrum for higher resolution

Input / Output Specifications¶

Input (oscillation): RheoData with x = \(\omega\) (rad/s), y = complex \(G^*\) or (N, 2) array \([G', G'']\)
Input (relaxation): RheoData with x = time (s), y = \(G(t)\) (Pa)
Output: RheoData with x = \(\tau\) (s), y = \(H(\tau)\) (Pa)

Metadata includes SpectrumResult with regularization_param, residual_norm, and method.

Usage¶

From Oscillatory Data¶

from rheojax.transforms import SpectrumInversion
from rheojax.core.data import RheoData
import numpy as np

# Frequency sweep data
omega = np.logspace(-2, 2, 50)
G_prime = 1e4 * omega**2 / (1 + omega**2)  # Maxwell model
G_double_prime = 1e4 * omega / (1 + omega**2)
G_star = G_prime + 1j * G_double_prime

data = RheoData(x=omega, y=G_star, metadata={'test_mode': 'oscillation'})

# Recover relaxation spectrum
inv = SpectrumInversion(method="tikhonov", n_tau=100)
spectrum_data, info = inv.transform(data)

tau = spectrum_data.x
H_tau = spectrum_data.y

# For a single Maxwell element, expect a peak at τ = 1 s
peak_tau = tau[np.argmax(H_tau)]
print(f"Peak relaxation time: {peak_tau:.3f} s")
print(f"Regularization λ: {info['spectrum_result'].regularization_param:.4g}")

From Relaxation Data¶

from rheojax.transforms import SpectrumInversion

# Stress relaxation data
t = np.logspace(-3, 2, 100)
G_t = 5e3 * np.exp(-t / 0.1) + 2e3 * np.exp(-t / 10.0)

relax_data = RheoData(x=t, y=G_t, metadata={'test_mode': 'relaxation'})

inv = SpectrumInversion(source="relaxation", method="max_entropy", n_tau=80)
spectrum_data, info = inv.transform(relax_data)

# Expect two peaks at τ ≈ 0.1 s and τ ≈ 10 s
import matplotlib.pyplot as plt
plt.semilogx(spectrum_data.x, spectrum_data.y)
plt.xlabel(r'$\tau$ (s)')
plt.ylabel(r'$H(\tau)$ (Pa)')
plt.title('Relaxation Spectrum')

Comparing Methods¶

from rheojax.transforms import SpectrumInversion

# Compare Tikhonov and MaxEnt on same data
tik = SpectrumInversion(method="tikhonov", n_tau=100)
mem = SpectrumInversion(method="max_entropy", n_tau=100)

spec_tik, info_tik = tik.transform(osc_data)
spec_mem, info_mem = mem.transform(osc_data)

print(f"Tikhonov residual: {info_tik['spectrum_result'].residual_norm:.4g}")
print(f"MaxEnt residual:   {info_mem['spectrum_result'].residual_norm:.4g}")

Validation and Quality Control¶

Diagnostic Checks¶

1. Residual norm: Measures data fidelity

Low residual + smooth spectrum → good inversion
Low residual + oscillatory spectrum → under-regularized (reduce \(\lambda\))
High residual + smooth spectrum → over-regularized (increase \(n_\tau\) or reduce \(\lambda\))

2. Back-calculation test:

# Verify: can H(τ) reproduce the original data?
from rheojax.transforms import PronyConversion
# Discretize H(τ) as Prony modes
G_i = H_tau * np.diff(np.log(tau), append=np.log(tau[-1]) - np.log(tau[-2]))
prony = PronyConversion(direction="time_to_freq")
# ... compare reconstructed G'(ω), G''(ω) with original

Common Failure Modes¶

1. Spurious oscillations:

Cause: Insufficient regularization (\(\lambda\) too small)
Fix: Increase regularization or switch to max_entropy

2. Over-smoothed spectrum:

Cause: Excessive regularization (\(\lambda\) too large)
Fix: Reduce regularization, increase n_tau

3. Truncated spectrum:

Cause: tau_range does not span the full relaxation window
Fix: Widen tau_range or let auto-detection handle it

4. Negative H(τ) values:

Cause: Tikhonov can produce negative values before clamping
Fix: Use max_entropy (inherently non-negative)

API References¶

Module: rheojax.transforms
Class: rheojax.transforms.SpectrumInversion

References¶

Honerkamp, J. & Weese, J. (1993). “A nonlinear regularization method for the calculation of relaxation spectra.” Rheol. Acta, 32, 65–73. DOI: 10.1007/BF00396678
Baumgaertel, M. & Winter, H.H. (1989). “Determination of discrete relaxation and retardation time spectra from dynamic mechanical data.” Rheol. Acta, 28, 511–519. DOI: 10.1007/BF01332922
Davies, A.R. & Anderssen, R.S. (1997). “Sampling localization in determining the relaxation spectrum.” J. Non-Newtonian Fluid Mech., 73, 163–179. DOI: 10.1016/S0377-0257(97)00056-6
Hansen, P.C. (1992). “Analysis of discrete ill-posed problems by means of the L-curve.” SIAM Review, 34, 561–580. DOI: 10.1137/1034115
Bryan, R.K. (1990). “Maximum entropy analysis of oversampled data problems.” Eur. Biophys. J., 18, 165–174. DOI: 10.1007/BF02427376
Tschoegl, N.W. (1989). The Phenomenological Theory of Linear Viscoelastic Behavior: An Introduction. Springer-Verlag. DOI: 10.1007/978-3-642-73602-5