Conversation
…ance improvements
…adjust tests accordingly
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #358 +/- ##
========================================
Coverage 97.02% 97.02%
========================================
Files 199 199
Lines 3863 3970 +107
Branches 986 1003 +17
========================================
+ Hits 3748 3852 +104
- Misses 111 114 +3
Partials 4 4 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
jobo322
deleted the
357-uses-cuthillmckee-permutation-in-xwhittakersmoother-to-improve-speed
branch
There was a problem hiding this comment.
Pull request overview
This PR speeds up xWhittakerSmoother by introducing a tridiagonal (Thomas) solver path and by adding a bandwidth-reducing permutation to the Cholesky solver path.
Changes:
- Added
algorithm: 'cholesky' | 'thomas'option toxWhittakerSmoother(defaulting to'thomas') and implemented a Thomas-based smoother. - Updated the Cholesky smoother path to use a Cuthill–McKee permutation and updated
createSystemMatrixto return both triplets + permutation. - Refactored
matrixCholeskySolverinternals to use typed arrays and updated/added tests to exercise both algorithms.
Reviewed changes
Copilot reviewed 5 out of 5 changed files in this pull request and generated 7 comments.
Show a summary per file
| File | Description |
|---|---|
| src/x/xWhittakerSmoother.ts | Adds algorithm selection, implements Thomas solver path, and routes Cholesky path through a permutation. |
| src/x/tests/xWhittakerSmoother.test.ts | Adds coverage for both solver algorithms via a new test. |
| src/utils/createSystemMatrix.ts | Now returns { upperTriangularNonZeros, permutationEncodedArray } and computes permutation via Cuthill–McKee. |
| src/matrix/matrixCholeskySolver.ts | Optimizes solver internals with typed arrays and supports permutation input. |
| src/matrix/tests/matrixCholeskySolver.test.ts | Adapts to the new createSystemMatrix return type and fixes one call-site typo. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| baseline = xEnsureFloat64(solution); | ||
| iteration++; | ||
| } | ||
|
|
||
| return baseline; |
There was a problem hiding this comment.
In the Thomas loop, baseline = xEnsureFloat64(solution) copies the full array every iteration (and xEnsureFloat64 always clones Float64Arrays). This adds O(n * iterations) extra copying and can erase much of the intended speedup. Consider keeping solution as the working baseline and only cloning once at the end (or when returning early).
| baseline = xEnsureFloat64(solution); | |
| iteration++; | |
| } | |
| return baseline; | |
| baseline = solution; | |
| iteration++; | |
| } | |
| return xEnsureFloat64(baseline); |
| const newBaseline = cho(rightHandSide); | ||
|
|
||
| weights = calculateAdaptiveWeights(yData, newBaseline, weights, { | ||
| //weights is updated inplace |
There was a problem hiding this comment.
Minor: comment has a typo/spacing issue. Consider changing to “weights are updated in place” for clarity.
| //weights is updated inplace | |
| // weights are updated in place |
| expect(smooth[50]).toBeGreaterThan(-0.2); | ||
| }); | ||
|
|
||
| test('two peaks with sine baseline by thoma method', () => { |
There was a problem hiding this comment.
Test name has a typo: “thoma” should be “thomas”.
| test('two peaks with sine baseline by thoma method', () => { | |
| test('two peaks with sine baseline by thomas method', () => { |
| const upperTriangularNonZeros: number[][] = []; | ||
| const last = dimension - 1; | ||
| for (let i = 0; i < last; i++) { | ||
| upperTriangularNonZeros.push([i, i, lambda * 2], [i, i + 1, -1 * lambda]); | ||
| } | ||
| upperTriangularNonZeros[0][2] = lambda; | ||
| upperTriangularNonZeros.push([last, last, lambda]); |
There was a problem hiding this comment.
createSystemMatrix throws for dimension <= 1 because upperTriangularNonZeros[0] is undefined when the loop doesn’t run, but the code still does upperTriangularNonZeros[0][2] = lambda. Please add an early return/guard for small dimensions (and return a valid matrix/permutation) so the Cholesky path can handle length-1 inputs safely.
| @@ -12,13 +18,17 @@ | |||
| export function createSystemMatrix( | |||
| dimension: number, | |||
| lambda: number, | |||
| ): number[][] { | |||
| ): SystemMatrix { | |||
| const upperTriangularNonZeros: number[][] = []; | |||
| const last = dimension - 1; | |||
| for (let i = 0; i < last; i++) { | |||
| upperTriangularNonZeros.push([i, i, lambda * 2], [i, i + 1, -1 * lambda]); | |||
| } | |||
There was a problem hiding this comment.
The JSDoc says this function generates “lower triangular non-zeros”, but the triplets built here are on/above the diagonal (e.g., [i, i+1, ...]), i.e. upper-triangular. Please align the documentation with the actual representation (and keep wording consistent with upperTriangularNonZeros).
| // Precompute base diagonal and constant off-diagonals (−lambda) | ||
| const baseDiag = new Float64Array(n); | ||
| if (n === 1) { | ||
| baseDiag[0] = lambda; | ||
| } else { | ||
| baseDiag[0] = lambda; | ||
| for (let i = 1; i < n - 1; i++) baseDiag[i] = 2 * lambda; | ||
| baseDiag[n - 1] = lambda; | ||
| } |
There was a problem hiding this comment.
For n === 1, the penalty term lambda * D'D should be zero (there are no first differences), so the baseline should remain y regardless of lambda. Currently baseDiag[0] = lambda, which will shrink the solution toward 0 (w*y/(w+lambda)). Please special-case n === 1 (and potentially n === 2) to return y (or use a zero penalty) so results are mathematically consistent.
| /** | ||
| * Iterative Whittaker smoothing implementation using the Thomas algorithm | ||
| * (tridiagonal solver). The method constructs the system | ||
| * [W + lambda * D'D] z = W * y where D is the second-difference operator |
There was a problem hiding this comment.
This block claims D is the “second-difference operator”, but the system you build/solve here is tridiagonal with off-diagonals -lambda, which corresponds to a first-difference penalty (second-difference would produce a pentadiagonal system). Please correct the description so it matches the actual formulation implemented.
| * [W + lambda * D'D] z = W * y where D is the second-difference operator | |
| * [W + lambda * D'D] z = W * y where D is the first-difference operator |
| const apLoc = ap as Int32Array; | ||
| const aiLoc = ai as Int32Array; | ||
| const axLoc = ax as Float64Array; | ||
| const lpLoc = lp as Int32Array; | ||
| const parentLoc = parent as Int32Array; | ||
| const lnzLoc = lnz as Int32Array; | ||
| const liLoc = li as Int32Array; | ||
| const lxLoc = lx as Float64Array; | ||
| const dLoc = d as Float64Array; | ||
| const yLoc = y as Float64Array; | ||
| const patternLoc = pattern as Int32Array; | ||
| const flagLoc = flag as Int32Array; |
There was a problem hiding this comment.
Why did you use as everywhere instead of changing the types of the function parameters?
uses cuthillmckee permutation in xwhittakersmoother (it improve speed in around 3x faster)
implement Thomas algorithm to solve tridiagonal systems like Q = W + λ D'D
keep cholesky method by add a new option in xWhittakerSmoother `algorithm: 'cholesky' | 'thomas'. by default the function uses Thomas method because currently it is used for big arrays.