Eight Failures That Never Errored: The Fitter Quietly Returned My Initial Guess

Writing an equivalent-circuit fitter for impedance spectra turned up eight problems. Not one threw an exception. Not one printed a warning. Every one of them produced output that looked entirely reasonable — and the worst of them had the fitting routine quietly hand back your initial guess while the interface reported success. This is not really about curve fitting. It is about how numerical code lies to you without erroring.

1. The expensive one: a sign error returns your initial guess

One step of Levenberg–Marquardt, with residuals defined as r = y − f(x) and Jacobian J = ∂r/∂x, descends along:

Δx = −(JᵀJ + λD)⁻¹ Jᵀr

I wrote:

x += (JᵀJ + λD)⁻¹ Jᵀr        // missing the minus sign

A missing minus sign means every step walks uphill. The inner LM loop responds by increasing the damping λ and retrying, rejecting any step whose residual is not better than the current one. Walking uphill, all twelve λ attempts get rejected, the loop exits, iter = 0, and the function returns the x it was handed on entry.

Which is the initial guess.

So the interface shows: a fitted curve (computed from the initial guess, roughly the right shape), a set of parameters (the ones you typed in), and a chi² value (finite, not absurd). Nothing anywhere indicates failure. With a decent starting guess the “result” even looks good.

The only way to catch it is to attack the fitter with synthetic data: forward-model a curve from a known parameter set, feed it in, and check whether those parameters come back out. If they do not, it is broken.

Every fitter, optimiser and solver needs a “synthetic data → known answer” self-test, run after every change to the numerical core. This is not a coverage question — code of this kind does not report its own failures, so you have to go and ask.

2. A parameter runs to 1e154 while the residual does not move

A ZARC element — a resistor in parallel with a constant phase element — has impedance:

Z = R / (1 + R·Q·(jω)ⁿ)

Push R toward infinity while holding the product R·Q fixed, and the expression degenerates into a pure CPE: Z → 1/(Q(jω)ⁿ).

Which means: along the curve where R·Q is constant, R and Q can run off to infinity together while the impedance barely changes. That is a flat valley, and the residual is the same everywhere along its floor. The optimiser has no reason to stop at any particular point, so it keeps sliding.

Without bounds, one of 40 noise realisations produced R1 ≈ 1e154 — with a perfectly normal chi². The fit “succeeded” and the parameters are garbage.

The remedy is not more iterations or a different optimiser (the valley is a property of the model, so every optimiser meets it), but two steps:

1. Give every parameter physical bounds (no resistance is 1e154 ohms);
2. When a parameter lands on a bound, mark it “unidentified” rather than printing the number.

The second step is the important one. Hitting a bound means the data does not constrain that parameter; the bound value is merely where the optimiser stopped and carries no physical meaning. Printing it dresses an arbitrary number up as a measurement.

3. A guard that cries wolf spends its own credibility

With bounds in place I added a reliability guard: a fit only counts as converged if the last iteration’s relative improvement is below 1e-12. That sounds rigorous.

It flagged 4 of 60 perfectly good fits as unreliable — fits whose reduced chi² sat in the same range as all the successful ones.

I had misread LM’s termination condition. The normal way LM finishes is by failing to find any step that reduces the residual — damping rises, nothing improves, the loop exits. At that point “the last relative improvement” may not exist at all, or may be the value from the previous accepted step. Using it as a convergence criterion misclassifies normal convergence as failure.

The cost is not those four fits. It is that a warning which fires wrongly teaches users to ignore it — and once that habit forms, the one time it fires correctly gets ignored too. A guard’s value rests entirely on rarely lying, and every false positive spends down that balance.

Under-report rather than over-report: the case you miss at least does not train the user to disregard you.

4. Weighting is not a setting, it is part of the result

Fitting an impedance spectrum requires choosing a residual weighting. Two common ones: unit weighting (every point equal) and modulus weighting (normalised by |Z|). Impedance magnitude can span several orders across a single spectrum, so this is not a detail.

Running the same synthetic data through 60 independent noise realisations and recording the recovery error on R0:

A 250x difference in systematic bias and a 160x difference in spread. Changing one dropdown gives two different conclusions from identical data.

So in this tool the weighting is an explicit control, and it is printed alongside the results. Hiding it behind a default hands the reader a number nobody can reproduce.

5. The remaining four: all in the display layer, all silent

A Nyquist plot must use equal axis scaling

Part of the point of reading a Nyquist plot is judging whether an arc is a proper semicircle or a depressed one. At a free aspect ratio, a depressed arc and a proper semicircle look identical — the exact thing you came to see is the thing you cannot see. Both axes must share one scale.

Do not generate ticks by accumulation

// bad: floating-point error accumulates
for (let v = lo; v <= hi; v += step) ticks.push(v);

// good: integer multiples
for (let k = 0; k <= n; k++) ticks.push(lo + k * step);

Accumulation turns a tick that should be exactly 0 into 1.7e-18. On its own that is harmless — but SI prefix formatting then renders it as 0.000001735 p, a string long enough to push the whole axis outside the viewBox. A 1e-18 numerical error presents as a broken chart.

Dimensionless quantities must not get SI prefixes

The CPE exponent n is dimensionless, typically around 0.86. A formatter applying SI prefixes indiscriminately displayed it as "859.6 m" — milli, as a prefix for 0.8596, is arithmetically correct and physically meaningless. The formatter needs to know which quantities carry units.

The test environment lies to you too

Twice during debugging I was certain a change had not taken effect, and both times it was the environment:

• headless Chrome reusing a --user-data-dir serves the old page out of its disk cache;
• a page URL without a cache-busting parameter serves the CDN edge copy — and note that the parameter name itself may be ignored by the cache rules. I used one that had no effect and misdiagnosed a whole round because of it.

Before verifying "did my change take effect," verify "am I looking at the new version."

6. What the eight have in common

Lined up together the pattern is unmistakable:

Not one of them raises an exception. This is how numerical code nearly always fails: it does not crash, it hands you a number. And how a number looks on screen is unrelated to whether it is right.

Only one thing defends against this — keep an input whose correct answer you already know, and run it after every change to the numerical core. The first four problems above were all found with synthetic data. None was found by reading the code.

References

• Marquardt (1963), An Algorithm for Least-Squares Estimation of Nonlinear Parameters
• Boukamp (1986), A nonlinear least squares fit procedure for analysis of immittance data
• impedance.py documentation: equivalent circuit fitting and weighting

Related: the EIS equivalent-circuit fitter (the tool described here), parameter fitting and identifiability (why "the fit converged" does not mean "the parameters are right"), and the HPPC pulse resistance analyser and dQ/dV incremental capacity analyser (the other two characterisation methods).