Relaxation spectra using nonlinear Tikhonov regularization with a Bayesian criterion

Abstract

Nonlinear Tikhonov regularization within a Bayesian framework is incorporated into a computer program called pyReSpect, which infers the continuous and discrete relaxation spectra from oscillatory shear experiments. It uses Bayesian inference to provide uncertainty estimates for the continuous spectrum h(τ) by propagating the uncertainty in the regularization parameter λ. The new algorithm is about 6–9 times faster than an older version of the program (ReSpect) in which the optimal λ was determined by the L-curve method. About half of the speedup arises from the Bayesian formulation by restricting the window of λ explored. The other half arises from the nonlinear formulation for which the spectrum is a weak function of λ, allowing us to use a coarse mesh for λ. The program is tested and validated on three examples: a synthetic spectrum, a H-polymer, and an elastomer with a nonzero terminal plateau.

Appendix: Linearity of the CRS with regularization parameter

What property of the problem leads to the general observation $\partial ^{2} H_{\lambda }/\partial (\log \lambda )^{2} \approx 0$? Since $H = \log h$, we intuitively understand that deviations in H are mild compared to deviations in h. Nevertheless, we can show using a toy example that the linearity of $H(\log \lambda )$ stems from the form of the cost function.

For simplicity, and without loss of generality, suppose that H is a variable rather than a function H(τ). Let, V (λ,H) = ρ²(H) + λη²(H) be a cost function that is linear in λ, but nonlinear in H. Furthermore, suppose

$$ \begin{array}{@{}rcl@{}} H^{+} & =& \min_{H} V(\lambda, H)\\ H^{*} & =& \min_{H} V(\lambda^{*}, H), \end{array} $$

where λ^∗ is any particular λ, not necessarily the optimal λ. Let R(H) = dρ²(H)/dH and N(H) = dη²(H)/dH, be the first derivatives. Optimality conditions imply

$$ \begin{array}{@{}rcl@{}} R(H^{+}) + \lambda N(H^{+}) & =& 0\\ R(H^{*}) + \lambda^{*} N(H^{*}) & =& 0. \end{array} $$

Let R⁺ = R(H⁺), N⁺ = N(H⁺), R^∗ = R(H^∗), and N^∗ = N(H^∗). For slowly varying functions, a Taylor series expansion implies

$$ \begin{array}{@{}rcl@{}} R^{+} & = R^{*} + (H^{+} - H^{*}) \left.\frac{dR}{dH}\right\vert_{H^{*}} + \mathcal{O}({\Delta} H^{2}) \end{array} $$

(40)

$$ \begin{array}{@{}rcl@{}} N^{+} & = N^{*} + (H^{+} - H^{*}) \left.\frac{dN}{dH}\right\vert_{H^{*}} + \mathcal{O}({\Delta} H^{2}) \end{array} $$

(41)

Ignoring quadratic and higher order terms, and setting c₁ = dR(H^∗)/dH and c₂ = dN(H^∗)/dH,

$$ \begin{array}{@{}rcl@{}} R^{+} + \lambda N^{+} & =& R^{*} + \lambda N^{*} + (H^{+} - H^{*}) (c_{1} + \lambda c_{2}) \\ 0 & =& (\lambda - \lambda^{*}) N^{*} + (H^{+} - H^{*}) (c_{1} + \lambda c_{2}). \end{array} $$

If c₂λ ≫ c₁, then (H⁺ − H^∗)/(λ − λ^∗) = ΔH/Δλ ≈−N^∗/(c₂λ). Thus,

$$ \frac{\Delta H}{\Delta \log \lambda} \approx -\frac{N^{*}}{c_{2}} = \text{ constant}. $$

(42)