diff --git a/apps/web/src/lib/document/presets.test.ts b/apps/web/src/lib/document/presets.test.ts index d1ec598..916ebaa 100644 --- a/apps/web/src/lib/document/presets.test.ts +++ b/apps/web/src/lib/document/presets.test.ts @@ -79,7 +79,8 @@ describe('graph presets', () => { 'sequences-summary', 'sequences-normalized-weights', 'sequences-first-seen', - 'sequences-linear-regression' + 'sequences-linear-regression', + 'sequences-logistic-regression' ]); expect( findGraphPreset('sequences-phyllotaxis')?.expressions.filter(({ source }) => @@ -98,6 +99,9 @@ describe('graph presets', () => { 'y = max(factors)', 'y = product(factors)' ]); + expect( + findGraphPreset('sequences-logistic-regression')?.expressions.map(({ source }) => source) + ).toContain('ys ~ l/(1 + exp(-k*(xs - x0)))'); expect( findGraphPreset('sequences-normalized-weights')?.expressions.some(({ source }) => source.includes('sum(weights)') diff --git a/apps/web/src/lib/document/presets.ts b/apps/web/src/lib/document/presets.ts index 3f48106..af9c66c 100644 --- a/apps/web/src/lib/document/presets.ts +++ b/apps/web/src/lib/document/presets.ts @@ -254,6 +254,19 @@ export const presetGroups: GraphPresetGroup[] = [ 'y = f(x)' ], [-5, 5, -7, 7] + ), + preset( + 'sequences-logistic-regression', + 'Logistic regression', + [ + 'xs = [-4, -3, -2, -1, 0, 1, 2, 3, 4]', + 'ys = [0.12, 0.28, 0.61, 1.39, 2.91, 4.32, 5.26, 5.72, 5.91]', + '(xs, ys)', + 'ys ~ l/(1 + exp(-k*(xs - x0)))', + 'f(u) = l/(1 + exp(-k*(u - x0)))', + 'y = f(x)' + ], + [-5, 5, 0, 6.5] ) ] }, diff --git a/apps/web/src/routes/docs/+page.svx b/apps/web/src/routes/docs/+page.svx index ca810f1..f7397d0 100644 --- a/apps/web/src/routes/docs/+page.svx +++ b/apps/web/src/routes/docs/+page.svx @@ -148,7 +148,7 @@ Improper integrals and symbolic antiderivatives are not currently supported. | slice | `values[2:4]` | second through fourth values | | indexed selection | `values[[1, 3]]` | first and third values | -### Linear regression +### Regression Use `~` to fit undefined scalar parameters in a model to finite sequence data: @@ -166,13 +166,26 @@ and root-mean-square error. Fitted parameters are document-level definitions, so later variables, functions, and plots can reuse them. An existing assignment, such as `b = 0`, fixes that value instead of fitting it. -The current regression solver accepts models that are linear in every fitted -parameter. This includes polynomial models such as -`ys ~ a*xs^2 + b*xs + c`: the powers of the data do not make the parameters -nonlinear. Models such as `ys ~ a*xs^b`, logistic curves, and sinusoidal fits -require nonlinear optimization and are not supported yet. Data sequences must -be nonempty and equally sized, and the supplied data must uniquely determine -every fitted parameter. +Models linear in every fitted parameter use a direct deterministic solver. This +includes polynomial models such as `ys ~ a*xs^2 + b*xs + c`: the powers of the +data do not make the parameters nonlinear. + +Genuinely nonlinear models use bounded nonlinear least squares. For example: + +```text +ys ~ a*exp(b*xs) + c +ys ~ a*xs^b +ys ~ l/(1 + exp(-k*(xs - x0))) +ys ~ a*sin(b*xs + c) + d +``` + +Nonlinear fitting is iterative and can have multiple local solutions, +particularly for sinusoidal models. The solver tries a deterministic set of +starting points and reports the best finite solution it finds; equivalent +parameterizations can therefore display different parameter values for the +same curve. A nonlinear relation may fit at most 8 parameters and 20,000 data +rows. All regression data must be nonempty and equally sized, and the supplied +data must independently determine every fitted parameter. Scalar arithmetic broadcasts over sequences. Equal-length sequences support elementwise `+` and `-`; `.*` is elementwise multiplication. A pair of diff --git a/crates/cartesium-core/src/scene/regression.rs b/crates/cartesium-core/src/scene/regression.rs index 1ae7e94..3af0333 100644 --- a/crates/cartesium-core/src/scene/regression.rs +++ b/crates/cartesium-core/src/scene/regression.rs @@ -9,6 +9,12 @@ use super::resolution::{ParsedExpressions, push_diagnostic}; use super::{GraphDocument2d, RegressionParameter, RegressionResult, SceneDiagnostic}; const MAX_REGRESSION_PARAMETERS: usize = 32; +const MAX_NONLINEAR_PARAMETERS: usize = 8; +const MAX_NONLINEAR_ROWS: usize = 20_000; +const MAX_NONLINEAR_EVALUATIONS: usize = 768; +const MAX_NONLINEAR_ITERATIONS_PER_START: usize = 80; + +type NonlinearCandidate = (Vec, Vec, f64); pub(super) fn raw_regression_parameter_names(parsed: &ParsedExpressions) -> BTreeSet { let assigned = parsed @@ -85,7 +91,7 @@ pub(super) fn resolve_regressions( if parameters.iter().any(|name| duplicates.contains(name)) { continue; } - match fit_linear_relation(&left, &right, ¶meters) { + match fit_relation(&left, &right, ¶meters) { Ok(result) => { for parameter in &result.parameters { variables.insert(parameter.name.clone(), Expr::Number(parameter.value)); @@ -130,13 +136,6 @@ fn resolve_relation( }) { return Err(format!("cannot fit reserved name '{name}'")); } - if !is_affine_in_parameters(&left, ¶meters) || !is_affine_in_parameters(&right, ¶meters) - { - return Err( - "nonlinear regression is not supported yet; fitted parameters must appear linearly" - .to_owned(), - ); - } Ok((left, right, parameters)) } @@ -145,6 +144,18 @@ struct FitResult { rmse: f64, } +fn fit_relation( + left: &Expr, + right: &Expr, + parameters: &BTreeSet, +) -> Result { + if is_affine_in_parameters(left, parameters) && is_affine_in_parameters(right, parameters) { + fit_linear_relation(left, right, parameters) + } else { + fit_nonlinear_relation(left, right, parameters) + } +} + fn fit_linear_relation( left: &Expr, right: &Expr, @@ -199,6 +210,353 @@ fn fit_linear_relation( }) } +fn fit_nonlinear_relation( + left: &Expr, + right: &Expr, + parameters: &BTreeSet, +) -> Result { + let parameter_names = parameters.iter().cloned().collect::>(); + if parameter_names.len() > MAX_NONLINEAR_PARAMETERS { + return Err(format!( + "nonlinear regression exceeds the maximum of {MAX_NONLINEAR_PARAMETERS} fitted parameters" + )); + } + + let mut evaluations = 0; + let mut best: Option = None; + let starts = nonlinear_starts(parameter_names.len()); + let screening_budget = MAX_NONLINEAR_EVALUATIONS * 45 / 100; + for (start_index, start) in starts.iter().cloned().enumerate() { + let remaining_starts = starts.len() - start_index; + let start_budget = (screening_budget.saturating_sub(evaluations) / remaining_starts) + .max(parameter_names.len() + 2); + let evaluation_limit = (evaluations + start_budget).min(screening_budget); + if evaluations + parameter_names.len() + 1 >= evaluation_limit { + break; + } + if let Some(candidate) = optimize_nonlinear( + left, + right, + ¶meter_names, + start, + &mut evaluations, + evaluation_limit, + )? && best + .as_ref() + .is_none_or(|(_, _, best_cost)| candidate.2 < *best_cost) + { + best = Some(candidate); + } + } + let Some((solution, fitted_residual, _)) = best else { + return Err("nonlinear regression could not find a finite starting point".to_owned()); + }; + let best_cost = squared_norm(&fitted_residual); + let refined = optimize_nonlinear( + left, + right, + ¶meter_names, + solution.clone(), + &mut evaluations, + MAX_NONLINEAR_EVALUATIONS - parameter_names.len(), + )?; + let (solution, fitted_residual, _) = match refined { + Some(refined) if refined.2 <= best_cost => refined, + _ => (solution, fitted_residual, best_cost), + }; + if fitted_residual.len() < parameter_names.len() { + return Err(format!( + "regression needs at least {} data rows for {} fitted parameters", + parameter_names.len(), + parameter_names.len() + )); + } + if fitted_residual.len() > MAX_NONLINEAR_ROWS { + return Err(format!( + "nonlinear regression exceeds the maximum of {MAX_NONLINEAR_ROWS} data rows" + )); + } + + // A local minimum is not a meaningful fit when the parameters are not + // independently observable there. Reuse the pivoted QR rank test rather + // than accepting an arbitrary point on a flat family of solutions. + let jacobian = finite_difference_jacobian( + left, + right, + ¶meter_names, + &solution, + &fitted_residual, + &mut evaluations, + )?; + solve_least_squares(jacobian, &vec![0.0; fitted_residual.len()])?; + + let rmse = (squared_norm(&fitted_residual) / fitted_residual.len() as f64).sqrt(); + Ok(FitResult { + parameters: parameter_names + .into_iter() + .zip(solution) + .map(|(name, value)| RegressionParameter { name, value }) + .collect(), + rmse, + }) +} + +fn nonlinear_starts(parameter_count: usize) -> Vec> { + let mut starts = vec![ + vec![0.0; parameter_count], + vec![1.0; parameter_count], + vec![-1.0; parameter_count], + (0..parameter_count) + .map(|index| if index % 2 == 0 { 1.0 } else { -1.0 }) + .collect(), + (0..parameter_count) + .map(|index| if index % 2 == 0 { -1.0 } else { 1.0 }) + .collect(), + ]; + for index in 0..parameter_count { + let mut scaled = vec![1.0; parameter_count]; + scaled[index] = 4.0; + starts.push(scaled); + } + starts +} + +fn optimize_nonlinear( + left: &Expr, + right: &Expr, + parameter_names: &[String], + mut values: Vec, + evaluations: &mut usize, + evaluation_limit: usize, +) -> Result, String> { + let Some(mut residual_values) = + try_residual(left, right, parameter_names, &values, evaluations)? + else { + return Ok(None); + }; + if residual_values.len() > MAX_NONLINEAR_ROWS { + return Err(format!( + "nonlinear regression exceeds the maximum of {MAX_NONLINEAR_ROWS} data rows" + )); + } + if residual_values.len() < parameter_names.len() { + return Err(format!( + "regression needs at least {} data rows for {} fitted parameters", + parameter_names.len(), + parameter_names.len() + )); + } + + let mut cost = squared_norm(&residual_values); + let mut damping = 1e-3; + for _ in 0..MAX_NONLINEAR_ITERATIONS_PER_START { + if *evaluations + parameter_names.len() + 1 >= evaluation_limit { + break; + } + let jacobian = finite_difference_jacobian( + left, + right, + parameter_names, + &values, + &residual_values, + evaluations, + )?; + let (normal, gradient) = normal_equations(&jacobian, &residual_values); + let gradient_scale = gradient.iter().map(|value| value.abs()).fold(0.0, f64::max); + if gradient_scale <= 1e-10 * cost.sqrt().max(1.0) { + break; + } + + let mut damped = normal; + for (index, row) in damped.iter_mut().enumerate() { + row[index] += damping * row[index].max(1.0); + } + let Some(step) = solve_square_system(damped, gradient.iter().map(|value| -value).collect()) + else { + damping *= 10.0; + continue; + }; + let trial = values + .iter() + .zip(&step) + .map(|(value, step)| value + step) + .collect::>(); + let Some(trial_residual) = try_residual(left, right, parameter_names, &trial, evaluations)? + else { + damping *= 10.0; + continue; + }; + if trial_residual.len() != residual_values.len() { + return Err("regression model changes the number of data rows".to_owned()); + } + let trial_cost = squared_norm(&trial_residual); + if trial_cost < cost { + let relative_step = step + .iter() + .zip(&values) + .map(|(step, value)| step.abs() / value.abs().max(1.0)) + .fold(0.0, f64::max); + values = trial; + residual_values = trial_residual; + let relative_improvement = (cost - trial_cost) / cost.max(1.0); + cost = trial_cost; + damping = (damping * 0.2).max(1e-12); + if relative_step < 1e-9 || relative_improvement < 1e-12 { + break; + } + } else { + damping = (damping * 10.0).min(1e16); + if damping >= 1e16 { + break; + } + } + } + Ok(Some((values, residual_values, cost))) +} + +fn finite_difference_jacobian( + left: &Expr, + right: &Expr, + parameter_names: &[String], + values: &[f64], + baseline: &[f64], + evaluations: &mut usize, +) -> Result>, String> { + let mut columns = Vec::with_capacity(values.len()); + for index in 0..values.len() { + if *evaluations >= MAX_NONLINEAR_EVALUATIONS { + return Err("nonlinear regression exceeded its evaluation budget".to_owned()); + } + let step = f64::EPSILON.cbrt() * values[index].abs().max(1.0); + let mut perturbed = values.to_vec(); + perturbed[index] += step; + let Some(perturbed_residual) = + try_residual(left, right, parameter_names, &perturbed, evaluations)? + else { + perturbed[index] = values[index] - step; + let Some(backward_residual) = + try_residual(left, right, parameter_names, &perturbed, evaluations)? + else { + return Err( + "nonlinear regression could not evaluate a parameter derivative".to_owned(), + ); + }; + if backward_residual.len() != baseline.len() { + return Err("regression model changes the number of data rows".to_owned()); + } + columns.push( + baseline + .iter() + .zip(backward_residual) + .map(|(baseline, backward)| (baseline - backward) / step) + .collect(), + ); + continue; + }; + if perturbed_residual.len() != baseline.len() { + return Err("regression model changes the number of data rows".to_owned()); + } + columns.push( + perturbed_residual + .into_iter() + .zip(baseline) + .map(|(perturbed, baseline)| (perturbed - baseline) / step) + .collect(), + ); + } + Ok(columns) +} + +fn try_residual( + left: &Expr, + right: &Expr, + parameter_names: &[String], + values: &[f64], + evaluations: &mut usize, +) -> Result>, String> { + *evaluations += 1; + let parameters = parameter_names + .iter() + .zip(values) + .map(|(name, value)| (name.as_str(), *value)) + .collect::>(); + match residual(left, right, ¶meters) { + Ok(values) if values.iter().all(|value| value.is_finite()) => Ok(Some(values)), + Ok(_) => Ok(None), + Err(message) + if message == "value must be finite" + || message == "sequence elements must be finite numbers" => + { + Ok(None) + } + Err(message) => Err(message), + } +} + +fn normal_equations(columns: &[Vec], residual: &[f64]) -> (Vec>, Vec) { + let count = columns.len(); + let mut normal = vec![vec![0.0; count]; count]; + let mut gradient = vec![0.0; count]; + for left in 0..count { + gradient[left] = dot(&columns[left], residual); + for right in 0..=left { + let value = dot(&columns[left], &columns[right]); + normal[left][right] = value; + normal[right][left] = value; + } + } + (normal, gradient) +} + +fn solve_square_system(mut matrix: Vec>, mut target: Vec) -> Option> { + let size = target.len(); + for pivot_index in 0..size { + let pivot = (pivot_index..size).max_by(|left, right| { + matrix[*left][pivot_index] + .abs() + .total_cmp(&matrix[*right][pivot_index].abs()) + })?; + matrix.swap(pivot_index, pivot); + target.swap(pivot_index, pivot); + let diagonal = matrix[pivot_index][pivot_index]; + if !diagonal.is_finite() || diagonal.abs() <= f64::EPSILON { + return None; + } + let pivot_row = matrix[pivot_index].clone(); + let pivot_target = target[pivot_index]; + for (row, row_target) in matrix + .iter_mut() + .skip(pivot_index + 1) + .zip(target.iter_mut().skip(pivot_index + 1)) + { + let factor = row[pivot_index] / diagonal; + for (value, pivot_value) in row + .iter_mut() + .skip(pivot_index) + .zip(pivot_row.iter().skip(pivot_index)) + { + *value -= factor * pivot_value; + } + *row_target -= factor * pivot_target; + } + } + let mut solution = vec![0.0; size]; + for row in (0..size).rev() { + let known = (row + 1..size) + .map(|column| matrix[row][column] * solution[column]) + .sum::(); + solution[row] = (target[row] - known) / matrix[row][row]; + } + solution + .iter() + .all(|value| value.is_finite()) + .then_some(solution) +} + +fn squared_norm(values: &[f64]) -> f64 { + dot(values, values) +} + fn residual( left: &Expr, right: &Expr, @@ -328,8 +686,8 @@ fn solve_least_squares(mut columns: Vec>, target: &[f64]) -> Result>, target: &[f64]) -> Result>(); + for (name, value) in expected { + assert!( + (parameters[name] - value).abs() < tolerance, + "{name}: expected {value}, got {}; all parameters: {parameters:?}", + parameters[name] + ); + } + } + + // Sine parameters have equivalent sign/phase representations, so test + // the recovered curve and identifiable magnitudes instead of one spelling. + let sinusoidal = build_scene_2d(&document(&[ + "xs = [-2, -1, 0, 1, 2, 3, 4]", + "ys = [-2.1346190772, -1.3105442181, 1.2735458558, 2.7937374665, 1.5887534296, -1.0245903523, -2.1904115221]", + "ys ~ a*sin(b*xs+c)+d", ])) .unwrap(); - assert_eq!(nonlinear.regressions.len(), 0); - assert_eq!(nonlinear.diagnostics.len(), 1); assert!( - nonlinear.diagnostics[0] - .message - .contains("nonlinear regression is not supported yet") + sinusoidal.diagnostics.is_empty(), + "{:?}", + sinusoidal.diagnostics ); + let regression = &sinusoidal.regressions[0]; + let parameters = regression + .parameters + .iter() + .map(|parameter| (parameter.name.as_str(), parameter.value)) + .collect::>(); + assert!((parameters["a"].abs() - 2.5).abs() < 2e-4); + assert!((parameters["b"].abs() - 1.1).abs() < 2e-4); + assert!((parameters["d"] - 0.3).abs() < 2e-4); + assert!(regression.rmse < 1e-8); +} +#[test] +fn regression_reports_underdetermined_models() { let rank_deficient = build_scene_2d(&document(&[ "xs = [1, 1, 1]", "ys = [2, 3, 4]", @@ -224,6 +284,34 @@ fn regression_reports_unsupported_and_underdetermined_models() { .message .contains("not uniquely determined") ); + + let nonlinear_rank_deficient = build_scene_2d(&document(&[ + "xs = [0, 0, 0]", + "ys = [2, 3, 4]", + "ys ~ a*exp(b*xs)", + ])) + .unwrap(); + assert_eq!(nonlinear_rank_deficient.regressions.len(), 0); + assert_eq!(nonlinear_rank_deficient.diagnostics.len(), 1); + assert!( + nonlinear_rank_deficient.diagnostics[0] + .message + .contains("not uniquely determined") + ); + + let no_finite_start = build_scene_2d(&document(&[ + "xs = [0, 1, 2]", + "ys = [1, 2, 3]", + "ys ~ ln(a*xs)", + ])) + .unwrap(); + assert_eq!(no_finite_start.regressions.len(), 0); + assert_eq!(no_finite_start.diagnostics.len(), 1); + assert!( + no_finite_start.diagnostics[0] + .message + .contains("could not find a finite starting point") + ); } #[test] diff --git a/docs/architecture.md b/docs/architecture.md index a4949d7..16a232f 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -116,14 +116,17 @@ remain scalar-only. Regression rows use `left ~ right` to fit undefined scalar names against finite sequence data. Existing assignments are fixed inputs; remaining names become document-level fitted parameters and can feed other assignments, functions, or -plots. The first solver accepts only expressions proven affine in every fitted -parameter. It evaluates the constant and unit-parameter bases through the -ordinary typed sequence evaluator, then solves the resulting least-squares -system with column-pivoted, reorthogonalized modified Gram-Schmidt QR. Empty or -mismatched data, rank-deficient systems, non-finite values, and nonlinear -parameter use produce localized diagnostics. Scene metadata reports the fitted -parameters and RMSE for the expression panel. Nonlinear optimization and table -input remain later slices. +plots. Models proven affine in every fitted parameter evaluate their constant +and unit-parameter bases through the ordinary typed sequence evaluator, then +solve the resulting least-squares system with column-pivoted, reorthogonalized +modified Gram-Schmidt QR. Models nonlinear in their fitted parameters use a +bounded multi-start Levenberg-Marquardt-style least-squares path with +finite-difference Jacobians. Parameter, row, iteration, and evaluation caps +bound scene-build work; a final pivoted-QR Jacobian check rejects locally +unidentifiable fits. The direct affine path remains preferred for its speed and +determinism. Empty or mismatched data, rank-deficient systems, and non-finite +models produce localized diagnostics. Scene metadata reports fitted parameters +and RMSE for the expression panel. Table input remains a later slice. The scalar elementary surface includes circular and inverse trigonometry (`sin`, `cos`, `tan`, `asin`, `acos`, `atan`, and `atan2`), reciprocal diff --git a/docs/project-status.md b/docs/project-status.md index 1733b11..81fd061 100644 --- a/docs/project-status.md +++ b/docs/project-status.md @@ -565,9 +565,10 @@ The current web application includes: - sequence-based least-squares regression through relations such as `ys ~ m*xs + b`, with undefined scalar parameters inferred from the model, deterministic QR fitting for models linear in every fitted parameter, - reusable fitted definitions, RMSE reporting, and localized diagnostics for - mismatched, underdetermined, rank-deficient, or nonlinear models; nonlinear - regression and table input remain future slices; + bounded multi-start nonlinear fitting for exponential, power, logistic, + sinusoidal, and other user-written models, reusable fitted definitions, RMSE + reporting, and localized diagnostics for mismatched, underdetermined, + rank-deficient, or non-finite models; table input remains a future slice; - Julia-style dense matrix literals (`[1 2; 3 4]`), one-based row/column indexing (`A[2, 1]`), matrix previews, scalar scaling, matrix-vector and matrix-matrix `*`, matrix transpose, and explicit vector `dot`/`norm`; @@ -749,7 +750,7 @@ as an independently drawn family of level lines. ### Broader calculator parity -Tables, movable points, labels, nonlinear regressions, folders/notes, accessibility +Tables, movable points, labels, folders/notes, accessibility polish, import/export formats, collaboration, and account features are not part of the current implemented milestone and do not yet have committed designs. Choose among them after the retained renderer is stable and the persistence