Repository navigation
Feat/root finding - #18
Open
backlundtransform wants to merge 8 commits into
Open
backlundtransform wants to merge 8 commits into
backlundtransform wants to merge 8 commits into
Conversation
The library had exactly one root finder, NumericExtensions.NewtonRaphson, and it could not be trusted. It ran a fixed 100 iterations with no convergence test and no early exit, divided by the derivative without checking it for zero - yielding a silent Inf or NaN where the tangent went flat - took no tolerance or iteration budget, and gave the caller no way to tell a converged answer from a diverged one. Each of those 100 iterations called DerivativeExtensions.Derivate, which allocates a Pascal matrix per call, so an easy root cost 100 matrix allocations and 200 function evaluations no matter how early it actually converged. Add Numerics/RootFinding/ with Bisection, Secant, Brent and a real Newton, each returning a RootResult carrying the estimate, whether it converged, the iteration count and the residual, so failure is reported rather than hidden. Brent is the default for a bracketed search and is exposed as the FindRoot extension. Newton stops on a vanishing or non-finite slope and on a non-finite iterate instead of iterating on NaN, and takes an analytic derivative through an overload. NewtonRaphson keeps its signature and delegates to the new Newton. Its default derivative reproduces Derivate at order 1 exactly - the same backward difference over the same 1e-5 step - so existing results are unchanged, but without the per-iteration allocation and with an early exit. Migrate the two hand-rolled Newton loops in the physics module, KeplerOrbit.TrueAnomaly and LagrangePoints.SolveCollinear, onto the shared implementation. Both were already correct; this removes the duplication. Also make NumericsTests.FindRoots assert the root it computes - it asserted Math.Abs(2) == 2 and never looked at the result. Suite: 1570 passing. The three failures in TimeserieValidationTests are unrelated and predate this branch: TestData/CsvTestDataGenerator formats decimals with the current culture, so on a Swedish machine 6.44 is written as 6,44 and collides with the CSV delimiter. They pass under DOTNET_SYSTEM_GLOBALIZATION_INVARIANT=1 and on CI.
AGENTS.md requires the section README to be updated when public API lands in that section, and the root-finding module added a namespace without one. Add a Root Finding section covering FindRoot, the four methods and when each applies, RootResult and EnsureConverged, and point the existing Newton-Raphson bullet at it.
Measured before the Del 1 migration, so the migration can be shown not to cost anything and the v4.4 performance work has something to beat. The numbers that bear on v4.3 directly: reusing an LU factorization is 57x faster than refactoring per solve at N=256 (119us against 6782us), which is the whole argument for the "factor once" rule; and Matrix.Inverse is the most expensive operation measured at 35.5ms, five times a full LU factor-and-solve, which is what the Kalman filters pay today by forming S.Inverse() explicitly. Four findings belong to v4.4 rather than v4.3 but were discovered here, so they are written down: every Matrix allocates a second n*n array for the identity field that almost nothing reads, doubling the memory cost of every matrix; new VectorN(double[]) clones its input, so SpMV allocates its result twice; the triple-loop multiply grows 70x for a 4x dimension against 64x for pure O(n^3), the rest being cache misses; and one MLP training epoch allocates 88.75 MB over 2048 samples. Runs use --inProcess: the Windows application control policy blocks the DLL copy BenchmarkDotNet places in its generated per-benchmark folder, so the default toolchain yields no results on this machine. Re-runs must match for the comparison to mean anything.
v4.1 built the decompositions but left every caller solving systems its own way. Five of those sites now go through LuDecomposition instead of their own Gaussian elimination: the GaussElimination extension itself, the RBF weight solve in MultivariateInterpolation, the polynomial regression solve in InferentialStatisticsExtensions, the dense not-a-knot fallback in CubicSplineInterpolation, and the global system in Assembler1D. The tridiagonal branch of the spline keeps the Thomas algorithm, which is O(n) against LU's O(n^3) and must not be replaced by it. GaussElimination was not merely duplicated, it was wrong. Its back substitution assigned the pivot A[i,i] to the solution component whenever the numerator came out exactly zero, so any system with a zero in its solution was answered incorrectly and silently: the diagonal system 2x=2, 3y=0, 4z=4 returned [1, 3, 1] instead of [1, 0, 1]. It now returns the right answer, and no longer modifies its arguments in place. The existing tests assert exact equality and still hold. Add LuDecomposition.SmallestPivotMagnitude so callers can detect a near-singular matrix rather than only an exactly singular one, and express IsSingular in terms of it. Without this the RBF solve would have lost its near-singularity guard and reported a bare "not invertible" in place of advice to change epsilon or kernel. Add characterisation tests for PanelMethod, which had neither tests nor callers. They assert invariants rather than physical correctness, because the one check against a closed-form solution does not pass: for a 200-panel cylinder at zero incidence the solver returns Cp ~= 1 on every panel, where potential flow gives a distribution over [-3, 1]. Whether the solver is wrong or a seamed cylinder is outside its domain is unresolved, so that test is left ignored with the question recorded, and PanelMethod itself is not migrated. Suite: 1580 passing, 1 ignored.
…sites Three migrations, continuing the work of routing callers onto the decompositions v4.1 built. The Kalman filters formed an explicit inverse of the innovation covariance to compute the gain. S is symmetric positive definite, so K = P*Ht*S^-1 is now solved as S*Kt = (P*Ht)t instead. The baseline measured Matrix.Inverse at 35.5ms for a 256x256 matrix against 6.8ms for a full LU factor-and-solve, and an explicit inverse can also break the symmetry the filter depends on. Add Matrix.SolveSymmetricPositiveDefinite, which uses Cholesky and falls back to LU if the matrix is not positive definite, since a recursively updated covariance can drift out of definiteness in floating point. It lives in LinearAlgebra rather than in StateEstimation because that is the section that owns the abstraction. All four inverse sites across KalmanFilter, ExtendedKalmanFilter and KalmanSmoother are converted. CoupledOscillators had its own private Jacobi eigensolver, duplicating EigenDecomposition. Its symmetrised dynamical matrix is symmetric by construction, so the shared decomposition takes its symmetric path and provides the same contract: ascending real eigenvalues, orthonormal eigenvectors as columns. PCA found its eigenpairs by power iteration with deflation. It now uses EigenDecomposition for both the covariance and the dual Gram path. This is an accuracy fix, not only de-duplication: two of the sixteen new PCA tests fail against the old implementation, with score columns correlated to 2.1e-5 and reconstruction off by 7.4e-8, which is what deflation error at a 1e-8 tolerance leaves behind. PowerIteration is deleted. MaxIterations, Tolerance and Seed no longer affect anything but are kept so existing calls and grid-search parameter dictionaries keep working; they are documented as inert. PCA had no tests at all before this. The new suite asserts properties of the decomposition rather than one implementation's output. Two planned items are dropped, with reasons recorded in the roadmap: PanelMethod is not migrated, and the legacy EigenValues, EigenVector and DominantEigenVector helpers do not delegate, because their published contract is rounded integer ratios that OdeSolver and existing tests depend on. Changing it is an API decision, not a refactor. Verification: the solution builds clean in Release on all three target frameworks, and the Kalman (6) and CoupledOscillators (45) tests passed before the Windows application control policy began blocking the test assembly continuously. The full suite has NOT been run against these three changes. Run it, or let CI do it, before merging.
…quations Every linear fitter built the design matrix and then immediately reduced it to XtX before solving. That squares the condition number of X, and a Vandermonde design - any polynomial fit - is poorly conditioned to begin with. The solve now happens on the design matrix itself. Measured on a degree-5 fit over [1, 2], by running this branch's new test against the previous implementation: the normal equations recover the coefficients to a worst-case error of 5.6e-5, QR reaches 7.2e-11. That is close to six decimal digits, and it is why the change is worth making rather than merely tidy. Add two primitives to FittingSolver: SolveLeastSquares, and SolveWeightedLeastSquares which scales each row by the square root of its weight so the weighted problem becomes an ordinary one. Both return (XtX)^-1 alongside the coefficients, recovered from the triangular factor as R^-1 R^-T, so the Gram matrix is never formed or inverted for standard errors either. LeastSquaresFitter, WeightedLeastSquaresFitter, RobustFitter (initial OLS, every IRLS step, and its covariance) and ParameterEstimation all move over, as does the covariance estimate in NonlinearLeastSquaresFitter. The Levenberg-Marquardt iteration itself does not move. Marquardt damping is applied to JtJ's diagonal, which has no design-matrix equivalent short of the augmented formulation [J; sqrt(lambda*D)] - a change to the algorithm, not to the solver, and so a separate decision. It is noted in the code, the roadmap and the section README. Rank deficiency is now detected by QrDecomposition.IsFullRank, which uses a relative threshold, rather than by an absolute 1e-14 pivot test. The exception type is unchanged, so NonlinearLeastSquaresFitter's catch still breaks its loop as before. Add FittingConditioningTests, which pin the accuracy down with a tolerance of 1e-9 - comfortably met by QR, far out of reach of the normal equations - so a revert to XtX fails loudly instead of quietly losing digits. Suite: 1602 passing, 1 ignored, verified in Release configuration.
Four items that had been on the list since v4.1 and never done. NaiveBayes.NumClasses threw NotImplementedException - the only such member left in the library, and the only IClassificationModel that did not provide it. It is now set during Fit as the highest label plus one, matching DecisionTree, RandomForest, KernelSVC, KNearestNeighbors and LinearSVC, so callers can use it to size label-indexed arrays. The model had no tests at all, so add a suite covering that, separation of two Gaussian blobs, classification of unseen points, the unfitted-model guard, and the hyperparameters-only Clone convention. Remove the csproj reference excluding NumericExtensions.cs~RF43b9e4.TMP. The file has not existed for a long time. Remove the unused SARSA._pendingNextAction field, which produced a CS0169 warning on every target framework. Library warnings drop from three to zero; the two that remain are unused locals in the test project and were not in scope. Consolidate four roadmaps that existed in both docs/ and docs/completed/. Three resolved straightforwardly, but the pairs were not simply duplicates and the choice mattered: the docs/ copy of AdvancedGameEngineRoadMap had become entirely unchecked while the completed/ copy was checked, and the completed/ copy is the one that matches the code - the Aerodynamics models are present, and AircraftState and FlightDynamicsEngine are absent only because Engines.Game was extracted into its own package. Multiphysics was the reverse: the docs/ copy recorded the 1-D FEM primitives as done, which they are, while the completed/ copy still called them deferred. That copy is replaced by the accurate one, with its Assembler1D line corrected to say LuDecomposition rather than Gaussian elimination, which this branch changed. Suite: 1608 passing, 1 ignored, verified in Release configuration.
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.