Skip to content

Commit 5d4a6f9

Browse files
soodokusoodoku
authored
The analytic Hessian was -2*grad/h; also fix cosine NW gradient and scale-dependent bandwidths (#10)
* Fix analytic Hessians, cosine NW weight derivative, and scale-dependent tolerance Three defects, each verified against central finite differences with a step-size convergence study before and after. 1. None of lscv, lscv_mv, loocv_mse or loocv_mse_mv returned a Hessian. In 1-D, d2/dh2 [h^-1 f(delta/h)] is h^-3 (2f + 4u f' + u^2 f''). The code used h^-2 (2f' + u f'') with a minus sign. Both summands of that term are odd in u, so it vanished exactly over the symmetric pair matrix, leaving hess == -2 * grad / h identically -- not a second derivative at all. For the Gaussian kernel the returned value was 5.879e-02 where the true second derivative is 2.563e-01. The multivariate versions had the analogous error: the bracket needs d(d+1) + 2(d+1)S + S^2 + sum_k U_k^2 (f''/f - (f'/f)^2), and the whole term belongs over h^(d+2), not h^(d+1). Same for the Nadaraya-Watson weight second derivative, w'' = (2K + 4uK' + u^2 K'')/h^3. It was wrong for gauss, biweight, triweight and cosine; epan was already right, and the Numba epan and unif kernels already used the correct general form. Finite differences now converge quadratically to the analytic value in every case: gauss lscv relative error 3.9e-04 -> 3.9e-06 -> 4.2e-08 as eps goes 1e-2 -> 1e-3 -> 1e-4. 2. The cosine NW weight first derivative was wrong. It read (pi/4)[(pi^2 u^2/4 - 1)cos - (pi u/2)sin]/h^2; the derivative of K(u)/h is (pi/4)[-cos + (pi u/2)sin]/h^2 -- a spurious quadratic term and the wrong sign on the sine. The LOOCV-MSE gradient was 1.833e-03 against a finite-difference value of -2.135e-03: wrong sign, so nw_bandwidth with kernel="cosine" walked the wrong way and returned Silverman's rule. 3. Bandwidths depended on the units of the data. The LSCV criterion is exactly scale-equivariant (c*LSCV(c*x, c*h) == LSCV(x, h) to 1.8e-16), so its minimiser must satisfy h(c*x) == c*h(x). Its gradient, however, scales like 1/c^2, and convergence was tested against an absolute tol=1e-5. For c >= 100 every kernel except unif reported convergence at the starting guess and returned Silverman's rule untouched -- 172% off for gauss. Tested on the dimensionless ratio |g|*h/|f| instead, which is invariant for both criteria. Equivariance now holds to 6e-12 or better across all kernels. Two consequences of the Hessian fix had to be handled: - With honest curvature the Newton step is dominated by the kinks that the C0 kernels put in the criterion, so the iterate crawled and the iteration cap bound long before the minimum: kde_bandwidth_mv with kernel="epan" fell to 0.47 against a true minimiser of 0.89. The line search now forward-tracks (doubles while the score keeps falling) as well as backtracking. Near a genuine quadratic minimum the first doubling is rejected, so local convergence is unaffected. - The Numba LOOCV scorers skipped any point whose weights had all underflowed, so a bandwidth small enough to underflow every weight scored exactly 0 -- a fake global minimum that the new line search could reach. They now agree with the NumPy reference (m = 0 there). Tests: each new test was confirmed to fail on the unfixed source. Suite: 77 passed before, 89 passed after. * Let a criterion that passes through zero converge _grad_converged scaled the gradient test by abs(f) and returned False outright when that scale was zero. An iterate sitting at f == 0 with a negligible gradient could therefore never converge, and burned the whole iteration budget with nothing left to gain. It now falls back to an absolute test at that single point. Also records, in _line_search's own docstring, that the accepted point is the last one that improved and so can lie past the nearest minimum -- the reason 1-D cosine settles on a local minimum 0.03% worse in objective than a dense grid. Refining that bracket was tried and rejected. Golden-sectioning it does find the minimum on a quadratic (1.2 against a true 1.0 becomes 1.004), but every variant needed an absolute floor or tolerance somewhere -- in the termination test, then in the h_floor clamp -- and each one broke h(c*x) == c*h(x) for unif at c=100 and cosine at c=1000. Reintroducing the scale dependence this branch exists to remove, in order to improve a line search that already returns a valid decrease, is a bad trade. The docstring now says what the invariant is for anyone who revisits it. Both raised by an independent review of the branch. Suite 89 -> 93. * Make optimizer scale invariant in both directions
1 parent fdab461 commit 5d4a6f9

7 files changed

Lines changed: 660 additions & 246 deletions

File tree

‎hbw/_numba_kde.py‎

Lines changed: 24 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -208,7 +208,7 @@ def _cosine_all(u: float) -> tuple[float, float, float, float, float, float]:
208208
sin_half = math.sin(_PI_2 * absu)
209209
k2 = _PI / 32.0 * (_PI * (2.0 - absu) * cos_half + 2.0 * sin_half)
210210
sign_u = 1.0 if u >= 0 else -1.0
211-
k2p = sign_u * (-_PI**3 / 64.0) * sin_half * (2.0 - absu)
211+
k2p = sign_u * (-(_PI**3) / 64.0) * sin_half * (2.0 - absu)
212212
k2pp = _PI**3 / 128.0 * (_PI * (absu - 2.0) * cos_half + 2.0 * sin_half)
213213
else:
214214
k2 = 0.0
@@ -235,12 +235,12 @@ def lscv_numba_gauss(x: NDArray, h: float) -> tuple[float, float, float]:
235235

236236
local_sums[tid, 0] += k2
237237
local_sums[tid, 1] += k2 + u * k2p
238-
local_sums[tid, 2] += 2.0 * k2p + u * k2pp
238+
local_sums[tid, 2] += 2.0 * k2 + 4.0 * u * k2p + u * u * k2pp
239239

240240
if i != j:
241241
local_sums[tid, 3] += k
242242
local_sums[tid, 4] += k + u * kp
243-
local_sums[tid, 5] += 2.0 * kp + u * kpp
243+
local_sums[tid, 5] += 2.0 * k + 4.0 * u * kp + u * u * kpp
244244

245245
sum_K2 = local_sums[:, 0].sum()
246246
sum_SF = local_sums[:, 1].sum()
@@ -256,12 +256,7 @@ def lscv_numba_gauss(x: NDArray, h: float) -> tuple[float, float, float]:
256256

257257
score = sum_K2 / (n2 * h) - 2.0 * sum_K / (nn1 * h)
258258
grad = -sum_SF / (n2 * h2) + 2.0 * sum_SK / (nn1 * h2)
259-
hess = (
260-
2.0 * sum_SF / (n2 * h3)
261-
- sum_SF2 / (n2 * h2)
262-
- 4.0 * sum_SK / (nn1 * h3)
263-
+ 2.0 * sum_SK2 / (nn1 * h2)
264-
)
259+
hess = sum_SF2 / (n2 * h3) - 2.0 * sum_SK2 / (nn1 * h3)
265260

266261
return score, grad, hess
267262

@@ -314,12 +309,12 @@ def lscv_numba_epan(x: NDArray, h: float) -> tuple[float, float, float]:
314309

315310
local_sums[tid, 0] += k2
316311
local_sums[tid, 1] += k2 + u * k2p
317-
local_sums[tid, 2] += 2.0 * k2p + u * k2pp
312+
local_sums[tid, 2] += 2.0 * k2 + 4.0 * u * k2p + u * u * k2pp
318313

319314
if i != j:
320315
local_sums[tid, 3] += k
321316
local_sums[tid, 4] += k + u * kp
322-
local_sums[tid, 5] += 2.0 * kp + u * kpp
317+
local_sums[tid, 5] += 2.0 * k + 4.0 * u * kp + u * u * kpp
323318

324319
sum_K2 = local_sums[:, 0].sum()
325320
sum_SF = local_sums[:, 1].sum()
@@ -335,12 +330,7 @@ def lscv_numba_epan(x: NDArray, h: float) -> tuple[float, float, float]:
335330

336331
score = sum_K2 / (n2 * h) - 2.0 * sum_K / (nn1 * h)
337332
grad = -sum_SF / (n2 * h2) + 2.0 * sum_SK / (nn1 * h2)
338-
hess = (
339-
2.0 * sum_SF / (n2 * h3)
340-
- sum_SF2 / (n2 * h2)
341-
- 4.0 * sum_SK / (nn1 * h3)
342-
+ 2.0 * sum_SK2 / (nn1 * h2)
343-
)
333+
hess = sum_SF2 / (n2 * h3) - 2.0 * sum_SK2 / (nn1 * h3)
344334

345335
return score, grad, hess
346336

@@ -393,12 +383,12 @@ def lscv_numba_unif(x: NDArray, h: float) -> tuple[float, float, float]:
393383

394384
local_sums[tid, 0] += k2
395385
local_sums[tid, 1] += k2 + u * k2p
396-
local_sums[tid, 2] += 2.0 * k2p + u * k2pp
386+
local_sums[tid, 2] += 2.0 * k2 + 4.0 * u * k2p + u * u * k2pp
397387

398388
if i != j:
399389
local_sums[tid, 3] += k
400390
local_sums[tid, 4] += k + u * kp
401-
local_sums[tid, 5] += 2.0 * kp + u * kpp
391+
local_sums[tid, 5] += 2.0 * k + 4.0 * u * kp + u * u * kpp
402392

403393
sum_K2 = local_sums[:, 0].sum()
404394
sum_SF = local_sums[:, 1].sum()
@@ -414,12 +404,7 @@ def lscv_numba_unif(x: NDArray, h: float) -> tuple[float, float, float]:
414404

415405
score = sum_K2 / (n2 * h) - 2.0 * sum_K / (nn1 * h)
416406
grad = -sum_SF / (n2 * h2) + 2.0 * sum_SK / (nn1 * h2)
417-
hess = (
418-
2.0 * sum_SF / (n2 * h3)
419-
- sum_SF2 / (n2 * h2)
420-
- 4.0 * sum_SK / (nn1 * h3)
421-
+ 2.0 * sum_SK2 / (nn1 * h2)
422-
)
407+
hess = sum_SF2 / (n2 * h3) - 2.0 * sum_SK2 / (nn1 * h3)
423408

424409
return score, grad, hess
425410

@@ -471,12 +456,12 @@ def lscv_numba_biweight(x: NDArray, h: float) -> tuple[float, float, float]:
471456

472457
local_sums[tid, 0] += k2
473458
local_sums[tid, 1] += k2 + u * k2p
474-
local_sums[tid, 2] += 2.0 * k2p + u * k2pp
459+
local_sums[tid, 2] += 2.0 * k2 + 4.0 * u * k2p + u * u * k2pp
475460

476461
if i != j:
477462
local_sums[tid, 3] += k
478463
local_sums[tid, 4] += k + u * kp
479-
local_sums[tid, 5] += 2.0 * kp + u * kpp
464+
local_sums[tid, 5] += 2.0 * k + 4.0 * u * kp + u * u * kpp
480465

481466
sum_K2 = local_sums[:, 0].sum()
482467
sum_SF = local_sums[:, 1].sum()
@@ -492,12 +477,7 @@ def lscv_numba_biweight(x: NDArray, h: float) -> tuple[float, float, float]:
492477

493478
score = sum_K2 / (n2 * h) - 2.0 * sum_K / (nn1 * h)
494479
grad = -sum_SF / (n2 * h2) + 2.0 * sum_SK / (nn1 * h2)
495-
hess = (
496-
2.0 * sum_SF / (n2 * h3)
497-
- sum_SF2 / (n2 * h2)
498-
- 4.0 * sum_SK / (nn1 * h3)
499-
+ 2.0 * sum_SK2 / (nn1 * h2)
500-
)
480+
hess = sum_SF2 / (n2 * h3) - 2.0 * sum_SK2 / (nn1 * h3)
501481

502482
return score, grad, hess
503483

@@ -559,12 +539,12 @@ def lscv_numba_triweight(x: NDArray, h: float) -> tuple[float, float, float]:
559539

560540
local_sums[tid, 0] += k2
561541
local_sums[tid, 1] += k2 + u * k2p
562-
local_sums[tid, 2] += 2.0 * k2p + u * k2pp
542+
local_sums[tid, 2] += 2.0 * k2 + 4.0 * u * k2p + u * u * k2pp
563543

564544
if i != j:
565545
local_sums[tid, 3] += k
566546
local_sums[tid, 4] += k + u * kp
567-
local_sums[tid, 5] += 2.0 * kp + u * kpp
547+
local_sums[tid, 5] += 2.0 * k + 4.0 * u * kp + u * u * kpp
568548

569549
sum_K2 = local_sums[:, 0].sum()
570550
sum_SF = local_sums[:, 1].sum()
@@ -580,12 +560,7 @@ def lscv_numba_triweight(x: NDArray, h: float) -> tuple[float, float, float]:
580560

581561
score = sum_K2 / (n2 * h) - 2.0 * sum_K / (nn1 * h)
582562
grad = -sum_SF / (n2 * h2) + 2.0 * sum_SK / (nn1 * h2)
583-
hess = (
584-
2.0 * sum_SF / (n2 * h3)
585-
- sum_SF2 / (n2 * h2)
586-
- 4.0 * sum_SK / (nn1 * h3)
587-
+ 2.0 * sum_SK2 / (nn1 * h2)
588-
)
563+
hess = sum_SF2 / (n2 * h3) - 2.0 * sum_SK2 / (nn1 * h3)
589564

590565
return score, grad, hess
591566

@@ -649,12 +624,12 @@ def lscv_numba_cosine(x: NDArray, h: float) -> tuple[float, float, float]:
649624

650625
local_sums[tid, 0] += k2
651626
local_sums[tid, 1] += k2 + u * k2p
652-
local_sums[tid, 2] += 2.0 * k2p + u * k2pp
627+
local_sums[tid, 2] += 2.0 * k2 + 4.0 * u * k2p + u * u * k2pp
653628

654629
if i != j:
655630
local_sums[tid, 3] += k
656631
local_sums[tid, 4] += k + u * kp
657-
local_sums[tid, 5] += 2.0 * kp + u * kpp
632+
local_sums[tid, 5] += 2.0 * k + 4.0 * u * kp + u * u * kpp
658633

659634
sum_K2 = local_sums[:, 0].sum()
660635
sum_SF = local_sums[:, 1].sum()
@@ -670,12 +645,7 @@ def lscv_numba_cosine(x: NDArray, h: float) -> tuple[float, float, float]:
670645

671646
score = sum_K2 / (n2 * h) - 2.0 * sum_K / (nn1 * h)
672647
grad = -sum_SF / (n2 * h2) + 2.0 * sum_SK / (nn1 * h2)
673-
hess = (
674-
2.0 * sum_SF / (n2 * h3)
675-
- sum_SF2 / (n2 * h2)
676-
- 4.0 * sum_SK / (nn1 * h3)
677-
+ 2.0 * sum_SK2 / (nn1 * h2)
678-
)
648+
hess = sum_SF2 / (n2 * h3) - 2.0 * sum_SK2 / (nn1 * h3)
679649

680650
return score, grad, hess
681651

@@ -751,24 +721,24 @@ def lscv_mv_numba_gauss(data: NDArray, h: float) -> tuple[float, float, float]:
751721
if k > 0:
752722
r = kp / k
753723
sum_ratio_k += u * r
754-
sum_d2_k += 2.0 * u * r + u * u * kpp / k
724+
sum_d2_k += u * u * (kpp / k - r * r)
755725

756726
if k2 > 0:
757727
r2 = k2p / k2
758728
sum_ratio_k2 += u * r2
759-
sum_d2_k2 += 2.0 * u * r2 + u * u * k2pp / k2
729+
sum_d2_k2 += u * u * (k2pp / k2 - r2 * r2)
760730

761731
local_sums[tid, 0] += prod_k2
762732
local_sums[tid, 1] += prod_k2 * (d + sum_ratio_k2)
763733
local_sums[tid, 2] += prod_k2 * (
764-
(d + 1) * d + 2.0 * (d + 1) * sum_ratio_k2 + sum_d2_k2
734+
(d + 1) * d + 2.0 * (d + 1) * sum_ratio_k2 + sum_ratio_k2 * sum_ratio_k2 + sum_d2_k2
765735
)
766736

767737
if i != j:
768738
local_sums[tid, 3] += prod_k
769739
local_sums[tid, 4] += prod_k * (d + sum_ratio_k)
770740
local_sums[tid, 5] += prod_k * (
771-
(d + 1) * d + 2.0 * (d + 1) * sum_ratio_k + sum_d2_k
741+
(d + 1) * d + 2.0 * (d + 1) * sum_ratio_k + sum_ratio_k * sum_ratio_k + sum_d2_k
772742
)
773743

774744
sum_K2 = local_sums[:, 0].sum()
@@ -786,12 +756,7 @@ def lscv_mv_numba_gauss(data: NDArray, h: float) -> tuple[float, float, float]:
786756

787757
score = sum_K2 / (n2 * hd) - 2.0 * sum_K / (nn1 * hd)
788758
grad = -sum_SF / (n2 * hd1) + 2.0 * sum_SK / (nn1 * hd1)
789-
hess = (
790-
(d + 1) * sum_SF / (n2 * hd2)
791-
- sum_SF2 / (n2 * hd1)
792-
- 2.0 * (d + 1) * sum_SK / (nn1 * hd2)
793-
+ 2.0 * sum_SK2 / (nn1 * hd1)
794-
)
759+
hess = sum_SF2 / (n2 * hd2) - 2.0 * sum_SK2 / (nn1 * hd2)
795760

796761
return score, grad, hess
797762

0 commit comments

Comments
 (0)