From 81146d1ec22cf4aec23f4f0bedb70fa0b900d1e2 Mon Sep 17 00:00:00 2001 From: tobiadefami Date: Thu, 10 Sep 2026 20:33:26 +0100 Subject: [PATCH 1/6] feat(HF-221): add LINEST regression function --- CHANGELOG.md | 4 + docs/guide/known-limitations.md | 8 + docs/guide/list-of-differences.md | 6 + src/error-message.ts | 2 + src/i18n/languages/csCZ.ts | 1 + src/i18n/languages/daDK.ts | 1 + src/i18n/languages/deDE.ts | 1 + src/i18n/languages/enGB.ts | 1 + src/i18n/languages/esES.ts | 1 + src/i18n/languages/fiFI.ts | 1 + src/i18n/languages/frFR.ts | 1 + src/i18n/languages/huHU.ts | 1 + src/i18n/languages/idID.ts | 1 + src/i18n/languages/itIT.ts | 1 + src/i18n/languages/nbNO.ts | 1 + src/i18n/languages/nlNL.ts | 1 + src/i18n/languages/plPL.ts | 1 + src/i18n/languages/ptPT.ts | 1 + src/i18n/languages/ruRU.ts | 1 + src/i18n/languages/svSE.ts | 1 + src/i18n/languages/trTR.ts | 1 + .../categories/statistical.ts | 12 ++ src/interpreter/plugin/RegressionPlugin.ts | 199 +++++++++++++++++ src/interpreter/plugin/index.ts | 1 + .../plugin/regression/LinearRegression.ts | 202 ++++++++++++++++++ 25 files changed, 451 insertions(+) create mode 100644 src/interpreter/plugin/RegressionPlugin.ts create mode 100644 src/interpreter/plugin/regression/LinearRegression.ts diff --git a/CHANGELOG.md b/CHANGELOG.md index aeea8dd511..5bad1b8ddb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), ## [Unreleased] +### Added + +- Added `LINEST` for simple and multiple linear regression, with optional intercept and regression statistics. The `stats` argument must be constant because result dimensions are determined before evaluation. + ### Fixed - Fixed the `AVERAGEIF` function returning a division-by-zero error when the calculated average was `0`. [#1733](https://github.com/handsontable/hyperformula/pull/1733) diff --git a/docs/guide/known-limitations.md b/docs/guide/known-limitations.md index 91b161d94a..3a6ceb1e4c 100644 --- a/docs/guide/known-limitations.md +++ b/docs/guide/known-limitations.md @@ -61,6 +61,14 @@ a circular reference. * Ordering (including mixed types, empty cells, and text collation) follows HyperFormula's own comparison rules, which honor the `caseSensitive` and `accentSensitive` configuration options. Numbers sort before text, and text before logical values. +### LINEST function + +* The `stats` argument controls whether the result has one or five rows. It must be omitted or supplied as a constant: `TRUE()`, `FALSE()`, a number, or the text `"TRUE"` or `"FALSE"`. Parentheses and numeric unary signs are supported. Cell references, named expressions, and calculated expressions such as `1=1` return `#VALUE!`, including when `LINEST` is nested inside another function. + +* Result width is determined from the predictor dimensions before evaluation. Changes to values within fixed input ranges and to a referenced `const` argument recalculate normally. Changes to dependency values do not resize the result automatically; this is the existing array-sizing limitation. + +* If an input expression returns a different predictor count than predicted, `LINEST` returns `#VALUE!`. For example, filtering predictor columns can change the required result width. + ### OFFSET function HyperFormula resolves the OFFSET function at parse time rather than during evaluation. The parser inspects the arguments and rewrites the expression into a plain cell reference or range. This keeps the dependency graph accurate but imposes several restrictions. diff --git a/docs/guide/list-of-differences.md b/docs/guide/list-of-differences.md index a109f82616..f53af063dc 100644 --- a/docs/guide/list-of-differences.md +++ b/docs/guide/list-of-differences.md @@ -34,6 +34,12 @@ See a full list of differences between HyperFormula, Microsoft Excel, and Google **Contents:** [[toc]] +## LINEST + +Unlike Excel, HyperFormula requires a constant `stats` argument because the output dimensions are determined before evaluation. For example, `=LINEST(A1:A10,B1:B10,TRUE(),D1)` returns `#VALUE!`; use `TRUE()` or `FALSE()` directly for `stats`. See [LINEST limitations](known-limitations.md#linest-function) for supported constants and input-sizing requirements. + +Numerical results can differ for nearly dependent predictors and nearly perfect fits. Coefficient standard errors are sensitive to conditioning, and the F statistic is sensitive to residuals close to machine precision. For an effectively perfect multiple regression, Excel and HyperFormula can return different large finite F values even when the coefficients agree. + ## General functionalities | Functionality | Examples | HyperFormula | Google Sheets | Microsoft Excel | diff --git a/src/error-message.ts b/src/error-message.ts index 4b80ee06fa..b431bdbc6c 100644 --- a/src/error-message.ts +++ b/src/error-message.ts @@ -7,6 +7,8 @@ * This is a class for detailed error messages across HyperFormula. */ export class ErrorMessage { + public static LinestStaticSize = 'LINEST requires input dimensions that determine a fixed result size.' + public static LinestStaticStats = 'LINEST requires a constant stats argument to determine its result size.' public static DistinctSigns = 'Distinct signs.' public static WrongArgNumber = 'Wrong number of arguments.' public static EmptyArg = 'Empty function argument.' diff --git a/src/i18n/languages/csCZ.ts b/src/i18n/languages/csCZ.ts index fea286083c..d3e0575403 100644 --- a/src/i18n/languages/csCZ.ts +++ b/src/i18n/languages/csCZ.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'JE.TEXT', LEFT: 'ZLEVA', LEN: 'DÉLKA', + LINEST: 'LINREGRESE', LN: 'LN', LOG10: 'LOG', LOG: 'LOGZ', diff --git a/src/i18n/languages/daDK.ts b/src/i18n/languages/daDK.ts index 9264551249..f3737cefae 100644 --- a/src/i18n/languages/daDK.ts +++ b/src/i18n/languages/daDK.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ER.TEKST', LEFT: 'VENSTRE', LEN: 'LÆNGDE', + LINEST: 'LINREGR', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/deDE.ts b/src/i18n/languages/deDE.ts index 4d3d39319b..0b1a180f0b 100644 --- a/src/i18n/languages/deDE.ts +++ b/src/i18n/languages/deDE.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ISTTEXT', LEFT: 'LINKS', LEN: 'LÄNGE', + LINEST: 'RGP', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/enGB.ts b/src/i18n/languages/enGB.ts index d271b732cd..3a7185fdad 100644 --- a/src/i18n/languages/enGB.ts +++ b/src/i18n/languages/enGB.ts @@ -144,6 +144,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ISTEXT', LEFT: 'LEFT', LEN: 'LEN', + LINEST: 'LINEST', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/esES.ts b/src/i18n/languages/esES.ts index 9c6b817b98..df2038ff36 100644 --- a/src/i18n/languages/esES.ts +++ b/src/i18n/languages/esES.ts @@ -143,6 +143,7 @@ export const dictionary: RawTranslationPackage = { ISTEXT: 'ESTEXTO', LEFT: 'IZQUIERDA', LEN: 'LARGO', + LINEST: 'ESTIMACION.LINEAL', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/fiFI.ts b/src/i18n/languages/fiFI.ts index 99ff05708f..87ab6be3cd 100644 --- a/src/i18n/languages/fiFI.ts +++ b/src/i18n/languages/fiFI.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ONTEKSTI', LEFT: 'VASEN', LEN: 'PITUUS', + LINEST: 'LINREGR', LN: 'LUONNLOG', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/frFR.ts b/src/i18n/languages/frFR.ts index f7b9be45f7..bc962747b3 100644 --- a/src/i18n/languages/frFR.ts +++ b/src/i18n/languages/frFR.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ESTTEXTE', LEFT: 'GAUCHE', LEN: 'NBCAR', + LINEST: 'DROITEREG', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/huHU.ts b/src/i18n/languages/huHU.ts index e40219553a..245709fa3b 100644 --- a/src/i18n/languages/huHU.ts +++ b/src/i18n/languages/huHU.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'SZÖVEG.E', LEFT: 'BAL', LEN: 'HOSSZ', + LINEST: 'LIN.ILL', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/idID.ts b/src/i18n/languages/idID.ts index 719cb11db6..532ee25e31 100644 --- a/src/i18n/languages/idID.ts +++ b/src/i18n/languages/idID.ts @@ -144,6 +144,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ADALAH.TEKS', LEFT: 'KIRI', LEN: 'PANJANG', + LINEST: 'LINEST', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/itIT.ts b/src/i18n/languages/itIT.ts index 1a494628f8..1a863ea46a 100644 --- a/src/i18n/languages/itIT.ts +++ b/src/i18n/languages/itIT.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'VAL.TESTO', LEFT: 'SINISTRA', LEN: 'LUNGHEZZA', + LINEST: 'REGR.LIN', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/nbNO.ts b/src/i18n/languages/nbNO.ts index ff15415a98..de4a1cf667 100644 --- a/src/i18n/languages/nbNO.ts +++ b/src/i18n/languages/nbNO.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ERTEKST', LEFT: 'VENSTRE', LEN: 'LENGDE', + LINEST: 'RETTLINJE', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/nlNL.ts b/src/i18n/languages/nlNL.ts index c491e2f44b..a7999b8913 100644 --- a/src/i18n/languages/nlNL.ts +++ b/src/i18n/languages/nlNL.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ISTEKST', LEFT: 'LINKS', LEN: 'PITUUS', + LINEST: 'LIJNSCH', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/plPL.ts b/src/i18n/languages/plPL.ts index 6fa4f016b4..7548158e88 100644 --- a/src/i18n/languages/plPL.ts +++ b/src/i18n/languages/plPL.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'CZY.TEKST', LEFT: 'LEWY', LEN: 'DŁ', + LINEST: 'REGLINP', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/ptPT.ts b/src/i18n/languages/ptPT.ts index 6c1e68c328..a74fba28b5 100644 --- a/src/i18n/languages/ptPT.ts +++ b/src/i18n/languages/ptPT.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ÉTEXTO', LEFT: 'ESQUERDA', LEN: 'NÚM.CARACT', + LINEST: 'PROJ.LIN', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/ruRU.ts b/src/i18n/languages/ruRU.ts index 6b35c1ce2f..e808973691 100644 --- a/src/i18n/languages/ruRU.ts +++ b/src/i18n/languages/ruRU.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ЕТЕКСТ', LEFT: 'ЛЕВСИМВ', LEN: 'ДЛСТР', + LINEST: 'ЛИНЕЙН', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/svSE.ts b/src/i18n/languages/svSE.ts index 5f87580087..7d869a89f1 100644 --- a/src/i18n/languages/svSE.ts +++ b/src/i18n/languages/svSE.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'ÄRTEXT', LEFT: 'VÄNSTER', LEN: 'LÄNGD', + LINEST: 'REGR', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/i18n/languages/trTR.ts b/src/i18n/languages/trTR.ts index 93c2176d37..991b50a31f 100644 --- a/src/i18n/languages/trTR.ts +++ b/src/i18n/languages/trTR.ts @@ -143,6 +143,7 @@ const dictionary: RawTranslationPackage = { ISTEXT: 'EMETİNSE', LEFT: 'SOL', LEN: 'UZUNLUK', + LINEST: 'DOT', LN: 'LN', LOG10: 'LOG10', LOG: 'LOG', diff --git a/src/interpreter/functionMetadata/categories/statistical.ts b/src/interpreter/functionMetadata/categories/statistical.ts index 1752c3eaf6..6a404977f5 100644 --- a/src/interpreter/functionMetadata/categories/statistical.ts +++ b/src/interpreter/functionMetadata/categories/statistical.ts @@ -304,6 +304,18 @@ export const STATISTICAL_DOCS: Record = { documentationUrl: 'https://hyperformula.handsontable.com/docs/guide/built-in-functions.html', examples: ['=LARGE(A1:A10, 1)', '=LARGE(A1:A10, 3)'], }, + LINEST: { + category: 'Statistical', + shortDescription: 'Returns linear regression coefficients and optional statistics.', + parameters: [ + {name: 'known_y', description: 'A numeric range of observed dependent values.'}, + {name: 'known_x', description: 'Optional numeric predictors. If omitted, uses sequential values starting at 1 with the shape of known_y.'}, + {name: 'const', description: 'Whether to fit an intercept. Defaults to TRUE; FALSE fits through zero.'}, + {name: 'stats', description: 'Whether to return five rows including regression statistics. Defaults to FALSE. Must be a constant; cell references and computed expressions are unsupported.'}, + ], + documentationUrl: 'https://hyperformula.handsontable.com/docs/guide/built-in-functions.html', + examples: ['=LINEST(A1:A10, B1:C10)', '=LINEST(A1:A10, B1:C10, TRUE(), TRUE())'], + }, 'LOGNORM.DIST': { category: 'Statistical', shortDescription: 'Returns density of lognormal distribution.', diff --git a/src/interpreter/plugin/RegressionPlugin.ts b/src/interpreter/plugin/RegressionPlugin.ts new file mode 100644 index 0000000000..d64c392281 --- /dev/null +++ b/src/interpreter/plugin/RegressionPlugin.ts @@ -0,0 +1,199 @@ +/** + * @license + * Copyright (c) 2025 Handsoncode. All rights reserved. + */ + +import {ArraySize} from '../../ArraySize' +import {CellError, ErrorType} from '../../Cell' +import {ErrorMessage} from '../../error-message' +import {Ast, AstNodeType, ProcedureAst} from '../../parser' +import {SimpleRangeValue} from '../../SimpleRangeValue' +import {coerceScalarToBoolean, coerceToRange} from '../ArithmeticHelper' +import {InterpreterState} from '../InterpreterState' +import {EmptyValue, getRawValue, InternalScalarValue, InterpreterValue, isExtendedNumber} from '../InterpreterValue' +import {FunctionArgumentType, FunctionPlugin, FunctionPluginTypecheck, ImplementedFunctions} from './FunctionPlugin' +import {fitLinearRegression, LinearRegressionResult} from './regression/LinearRegression' + +/** The orientation and predictor count shared by size prediction and runtime validation. */ +interface RegressionShape { + predictors: number, + orientation: 'paired' | 'columns' | 'rows', +} + +/** Classifies the original dimensions before any values are flattened. */ +function regressionShape(y: ArraySize, x?: ArraySize): RegressionShape | undefined { + if (x === undefined || (y.width === x.width && y.height === x.height)) { + return {predictors: 1, orientation: 'paired'} + } + if (y.width === 1 && y.height === x.height) { + return {predictors: x.width, orientation: 'columns'} + } + if (y.height === 1 && y.width === x.width) { + return {predictors: x.height, orientation: 'rows'} + } + return undefined +} + +/** LINEST rejects empty strings as Boolean options, including formula-generated strings. */ +function regressionBoolean(value: InternalScalarValue): boolean | CellError | undefined { + return value === '' ? undefined : coerceScalarToBoolean(value) +} + +/** Recognizes a numeric literal with optional parentheses and unary signs. */ +function isNumericConstant(ast: Ast): boolean { + if (ast.type === AstNodeType.PARENTHESIS) { + return isNumericConstant(ast.expression) + } + if (ast.type === AstNodeType.PLUS_UNARY_OP || ast.type === AstNodeType.MINUS_UNARY_OP) { + return isNumericConstant(ast.value) + } + return ast.type === AstNodeType.NUMBER +} + +/** Resolves scalar constants without evaluating dependencies during array-size prediction. */ +function staticBoolean(ast: Ast | undefined): boolean | undefined { + if (ast === undefined || ast.type === AstNodeType.EMPTY) { + return false + } + if (ast.type === AstNodeType.PARENTHESIS) { + return staticBoolean(ast.expression) + } + if (ast.type === AstNodeType.NUMBER || ast.type === AstNodeType.STRING) { + const value = regressionBoolean(ast.value) + return typeof value === 'boolean' ? value : undefined + } + if (ast.type === AstNodeType.PLUS_UNARY_OP || ast.type === AstNodeType.MINUS_UNARY_OP) { + return isNumericConstant(ast.value) ? staticBoolean(ast.value) : undefined + } + if (ast.type === AstNodeType.FUNCTION_CALL && ast.args.length === 0) { + if (ast.procedureName === 'TRUE') { + return true + } + if (ast.procedureName === 'FALSE') { + return false + } + } + return undefined +} + +/** Converts numerical failures to spreadsheet errors, including inside an array result. */ +function finiteValue(value: number): number | CellError { + return Number.isFinite(value) ? value : new CellError(ErrorType.NUM, ErrorMessage.NaN) +} + +/** Formats the five-row result while preserving statistics errors independently of coefficients. */ +function regressionOutput(fit: LinearRegressionResult, statistics: boolean, fitIntercept: boolean): SimpleRangeValue { + const result: InternalScalarValue[][] = [[...fit.coefficients].reverse().concat(fit.intercept).map(finiteValue)] + if (!statistics) { + return SimpleRangeValue.onlyValues(result) + } + const regressionSumSquares = fit.totalSumSquares - fit.residualSumSquares + const variance = fit.degreesOfFreedom === 0 ? 0 : fit.residualSumSquares / fit.degreesOfFreedom + const rSquared = fit.totalSumSquares === 0 ? 1 : regressionSumSquares / fit.totalSumSquares + const f = finiteValue((regressionSumSquares / fit.retainedPredictorCount) / variance) + const interceptError = fitIntercept ? finiteValue(fit.interceptError) : new CellError(ErrorType.NA) + result.push([...fit.standardErrors].reverse().map(finiteValue).concat(interceptError)) + result.push([finiteValue(rSquared), finiteValue(Math.sqrt(variance))]) + result.push([f, fit.degreesOfFreedom]) + result.push([finiteValue(regressionSumSquares), finiteValue(fit.residualSumSquares)]) + for (const row of result) { + while (row.length < fit.coefficients.length + 1) { + row.push(new CellError(ErrorType.NA)) + } + } + return SimpleRangeValue.onlyValues(result) +} + +/** Implements linear regression with statically predictable result dimensions. */ +export class RegressionPlugin extends FunctionPlugin implements FunctionPluginTypecheck { + public static implementedFunctions: ImplementedFunctions = { + 'LINEST': { + method: 'linest', + sizeOfResultArrayMethod: 'linestArraySize', + vectorizationForbidden: true, + enableArrayArithmeticForArguments: true, + parameters: [ + {argumentType: FunctionArgumentType.RANGE}, + {argumentType: FunctionArgumentType.ANY, defaultValue: EmptyValue}, + {argumentType: FunctionArgumentType.SCALAR, defaultValue: true, emptyAsDefault: true}, + {argumentType: FunctionArgumentType.SCALAR, defaultValue: false, emptyAsDefault: true}, + ], + }, + } + + /** + * Evaluates LINEST(known_y, [known_x], [const], [stats]). + * Returns coefficients in reverse predictor order, followed by the intercept. + * With stats enabled, adds standard errors; R-squared and residual standard error; + * F and residual degrees of freedom; regression and residual sums of squares. + * Unused statistics cells and the intercept standard error for const=FALSE are #N/A. + * + * A dynamic stats option is rejected here as well as during size prediction so that + * nested consumers cannot accidentally bypass the fixed-size contract. + */ + public linest(ast: ProcedureAst, state: InterpreterState): InterpreterValue { + const generatedX = ast.args.length < 2 || ast.args[1].type === AstNodeType.EMPTY + return this.runFunction(ast.args, state, this.metadata('LINEST'), + (knownY: SimpleRangeValue, knownX: InterpreterValue, constArg: InternalScalarValue, statsArg: InternalScalarValue) => { + const fitIntercept = regressionBoolean(constArg) + const statistics = regressionBoolean(statsArg) + if (fitIntercept instanceof CellError) { + return fitIntercept + } + if (statistics instanceof CellError) { + return statistics + } + if (fitIntercept === undefined || statistics === undefined) { + return new CellError(ErrorType.VALUE, ErrorMessage.WrongType) + } + if (staticBoolean(ast.args[3]) === undefined) { + return new CellError(ErrorType.VALUE, ErrorMessage.LinestStaticStats) + } + const x = generatedX ? undefined : coerceToRange(knownX) + const shape = regressionShape(knownY.size, x?.size) + if (shape === undefined) { + return new CellError(ErrorType.REF, ErrorMessage.ArrayDimensions) + } + if (shape.predictors + 1 > this.config.maxColumns) { + return new CellError(ErrorType.VALUE, ErrorMessage.ValueLarge) + } + if (this.linestArraySize(ast, state).width !== shape.predictors + 1) { + return new CellError(ErrorType.VALUE, ErrorMessage.LinestStaticSize) + } + const yValues = Array.from(knownY.valuesFromTopLeftCorner()) + const xValues = x === undefined ? yValues.map((_, i) => i + 1) : Array.from(x.valuesFromTopLeftCorner()) + if (!yValues.every(isExtendedNumber) || !xValues.every(isExtendedNumber)) { + return new CellError(ErrorType.VALUE, ErrorMessage.NumberRange) + } + const observations = yValues.map(getRawValue) as number[] + const numericX = xValues.map(getRawValue) as number[] + const predictors = observations.map((_, i) => { + if (shape.orientation === 'paired') { + return [numericX[i]] + } + return Array.from({length: shape.predictors}, (_, j) => numericX[shape.orientation === 'columns' ? i * shape.predictors + j : j * observations.length + i]) + }) + const fit = fitLinearRegression(predictors, observations, fitIntercept, statistics) + return regressionOutput(fit, statistics, fitIntercept) + }) + } + + /** Predicts output width from predictor geometry and height from a constant stats option. */ + public linestArraySize(ast: ProcedureAst, state: InterpreterState): ArraySize { + if (ast.args.length < 1 || ast.args.length > 4) { + return ArraySize.error() + } + const statistics = staticBoolean(ast.args[3]) + if (statistics === undefined) { + return ArraySize.error() + } + const arrayState = new InterpreterState(state.formulaAddress, true) + const y = this.arraySizeForAst(ast.args[0], arrayState) + const x = ast.args.length < 2 || ast.args[1].type === AstNodeType.EMPTY ? undefined : this.arraySizeForAst(ast.args[1], arrayState) + const shape = regressionShape(y, x) + if (shape === undefined || shape.predictors + 1 > this.config.maxColumns) { + return ArraySize.error() + } + return new ArraySize(shape.predictors + 1, statistics ? 5 : 1) + } +} diff --git a/src/interpreter/plugin/index.ts b/src/interpreter/plugin/index.ts index 2a3005375b..505776ea2b 100644 --- a/src/interpreter/plugin/index.ts +++ b/src/interpreter/plugin/index.ts @@ -51,3 +51,4 @@ export {StatisticalPlugin} from './StatisticalPlugin' export {MathPlugin} from './MathPlugin' export {ComplexPlugin} from './ComplexPlugin' export {StatisticalAggregationPlugin} from './StatisticalAggregationPlugin' +export {RegressionPlugin} from './RegressionPlugin' diff --git a/src/interpreter/plugin/regression/LinearRegression.ts b/src/interpreter/plugin/regression/LinearRegression.ts new file mode 100644 index 0000000000..d3a286db0b --- /dev/null +++ b/src/interpreter/plugin/regression/LinearRegression.ts @@ -0,0 +1,202 @@ +/** + * @license + * Copyright (c) 2025 Handsoncode. All rights reserved. + */ + +/** Numerical fit and uncertainties, before spreadsheet output formatting. */ +export interface LinearRegressionResult { + coefficients: number[], + intercept: number, + standardErrors: number[], + interceptError: number, + residualSumSquares: number, + totalSumSquares: number, + degreesOfFreedom: number, + retainedPredictorCount: number, +} + +/** A predictor column with its original position and pre-factorization norm. */ +interface RegressionColumn { + values: number[], + index: number, + originalNorm: number, +} + +/** + * Rank threshold relative to a predictor's original centered norm. + * Excel Online retains a 1e-5 perturbation of a duplicated predictor and removes 1e-6. + * Keep the threshold separate from roundoff handling for the fitted statistics. + */ +const RANK_TOLERANCE = 1e-6 + +/** Computes a Euclidean norm without squaring large inputs or spreading an unbounded range. */ +function norm(values: number[], start = 0): number { + let result = 0 + for (let i = start; i < values.length; i++) { + result = Math.hypot(result, values[i]) + } + return result +} + +/** Computes the mean relative to the first value to preserve small changes on large offsets. */ +function mean(values: number[]): number { + const origin = values[0] + let sum = 0 + for (const value of values) { + sum += (value - origin) / values.length + } + return origin + sum +} + +/** Applies a Householder reflection to the trailing part of a column. */ +function reflect(values: number[], vector: number[], start: number, factor: number): void { + let product = 0 + for (let i = start; i < values.length; i++) { + product += vector[i - start] * values[i] + } + product *= factor + for (let i = start; i < values.length; i++) { + values[i] -= product * vector[i - start] + } +} + +/** Solves the retained triangular system, treating an exactly zero predictor as a zero coefficient. */ +function backSubstitute(columns: RegressionColumn[], count: number, right: number[]): number[] { + const solution = Array(count).fill(0) + for (let i = count - 1; i >= 0; i--) { + if (columns[i].values[i] === 0) { + continue + } + let value = right[i] + for (let j = i + 1; j < count; j++) { + value -= columns[j].values[i] * solution[j] + } + solution[i] = value / columns[i].values[i] + } + return solution +} + +/** + * Fits observations using pivoted Householder QR, with a fixed leading intercept column. + * Predictor centering preserves accuracy for large offsets. The factorization is column-major, + * uses O(n*k + k*k) storage, and costs O(n*k*k) for n >= k; Q is never materialized. + * + * Exactly zero predictors intentionally consume an available QR slot. Excel Online returns a + * zero coefficient but still excludes that transformed response entry from residual statistics. + * A dependent nonzero predictor is instead removed using its relative norm, increasing df. + * + * @param predictors - observation rows, with predictor columns in their original order + * @param observations - numeric response for each observation + * @param fitIntercept - whether to include a constant term + * @param statistics - whether to compute coefficient uncertainties + */ +export function fitLinearRegression(predictors: number[][], observations: number[], fitIntercept: boolean, statistics: boolean): LinearRegressionResult { + const n = observations.length + const k = predictors[0].length + const offset = fitIntercept ? 1 : 0 + const means = Array(k).fill(0) + const columns: RegressionColumn[] = [] + if (fitIntercept) { + columns.push({values: Array(n).fill(1), index: -1, originalNorm: Math.sqrt(n)}) + } + for (let j = 0; j < k; j++) { + const values = predictors.map(row => row[j]) + means[j] = fitIntercept ? mean(values) : 0 + const centered = values.map(value => value - means[j]) + columns.push({values: centered, index: j, originalNorm: norm(centered)}) + } + + // Treat the intercept as the last input column moved to the front by a swap. + // Retain this order when pivot norms tie, so duplicated predictors agree with Excel. + if (fitIntercept && k > 1) { + const firstPredictor = columns.splice(1, 1)[0] + columns.push(firstPredictor) + } + // Remove a common response offset before reflection to preserve small variations. + const responseOrigin = fitIntercept ? observations[0] : 0 + const transformedY = observations.map(value => value - responseOrigin) + let remaining = columns.length + let count = 0 + while (count < Math.min(n, remaining)) { + if (count >= offset) { + let selected = count + let largest = -1 + for (let j = count; j < remaining; j++) { + const magnitude = norm(columns[j].values, count) + if (magnitude > largest) { + largest = magnitude + selected = j + } + } + [columns[count], columns[selected]] = [columns[selected], columns[count]] + } + const column = columns[count] + const magnitude = norm(column.values, count) + if (count >= offset && magnitude < RANK_TOLERANCE * column.originalNorm) { + remaining-- + columns[count] = columns[remaining] + columns[remaining] = column + continue + } + if (magnitude !== 0) { + const sign = column.values[count] >= 0 ? 1 : -1 + const vector = column.values.slice(count).map(value => value / magnitude) + vector[0] += sign + const factor = 1 / (1 + Math.abs(column.values[count]) / magnitude) + for (let j = count + 1; j < remaining; j++) { + reflect(columns[j].values, vector, count, factor) + } + reflect(transformedY, vector, count, factor) + column.values[count] = -sign * magnitude + column.values.fill(0, count + 1) + } + count++ + } + + const solution = backSubstitute(columns, count, transformedY) + const coefficients = Array(k).fill(0) + for (let j = offset; j < count; j++) { + coefficients[columns[j].index] = solution[j] + } + const intercept = fitIntercept ? mean(observations) - coefficients.reduce((sum, value, j) => sum + value * means[j], 0) : 0 + const yMean = fitIntercept ? mean(observations) : 0 + const centeredY = observations.map(value => value - yMean) + const totalSumSquares = norm(centeredY) ** 2 + let residualSumSquares = norm(transformedY, count) ** 2 + // Simple exact fits have zero residual statistics in Excel. Discard only + // reflection roundoff, whose squared scale is O((n * machine epsilon)^2). + const simpleFitRoundoff = k === 1 && residualSumSquares <= totalSumSquares * (n * Number.EPSILON) ** 2 + if (totalSumSquares === 0 || count === n || simpleFitRoundoff) { + residualSumSquares = 0 + } + const degreesOfFreedom = n - count + const variance = degreesOfFreedom === 0 ? 0 : residualSumSquares / degreesOfFreedom + const standardErrors = Array(k).fill(0) + let interceptVariance = 0 + if (statistics && variance !== 0) { + // Each solve produces one column of R^-1. Accumulate only covariance diagonals + // and the intercept's quadratic form, avoiding a full covariance matrix. + for (let j = 0; j < count; j++) { + const right = Array(count).fill(0) + right[j] = 1 + const inverseColumn = backSubstitute(columns, count, right) + let interceptWeight = fitIntercept ? inverseColumn[0] : 0 + for (let i = offset; i < count; i++) { + const original = columns[i].index + standardErrors[original] += variance * inverseColumn[i] ** 2 + interceptWeight -= means[original] * inverseColumn[i] + } + interceptVariance += variance * interceptWeight ** 2 + } + } + return { + coefficients, + intercept, + standardErrors: standardErrors.map(Math.sqrt), + interceptError: Math.sqrt(interceptVariance), + residualSumSquares, + totalSumSquares, + degreesOfFreedom, + retainedPredictorCount: count - offset, + } +} From 64b6393a8b88a0029be729a6554173243cca1ff7 Mon Sep 17 00:00:00 2001 From: tobiadefami Date: Thu, 10 Sep 2026 20:37:24 +0100 Subject: [PATCH 2/6] docs(HF-221): link LINEST changelog to PR 1769 --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5bad1b8ddb..0b5938d44a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,7 +9,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), ### Added -- Added `LINEST` for simple and multiple linear regression, with optional intercept and regression statistics. The `stats` argument must be constant because result dimensions are determined before evaluation. +- Added `LINEST` for simple and multiple linear regression, with optional intercept and regression statistics. The `stats` argument must be constant because result dimensions are determined before evaluation. [#1769](https://github.com/handsontable/hyperformula/pull/1769) ### Fixed From 65cb77397e606909f8245a7c3539a3e67f183645 Mon Sep 17 00:00:00 2001 From: tobiadefami Date: Thu, 10 Sep 2026 20:53:44 +0100 Subject: [PATCH 3/6] fix(HF-221): remove redundant LINEST type assertions --- src/interpreter/plugin/RegressionPlugin.ts | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/interpreter/plugin/RegressionPlugin.ts b/src/interpreter/plugin/RegressionPlugin.ts index d64c392281..7810eb2e7c 100644 --- a/src/interpreter/plugin/RegressionPlugin.ts +++ b/src/interpreter/plugin/RegressionPlugin.ts @@ -165,8 +165,8 @@ export class RegressionPlugin extends FunctionPlugin implements FunctionPluginTy if (!yValues.every(isExtendedNumber) || !xValues.every(isExtendedNumber)) { return new CellError(ErrorType.VALUE, ErrorMessage.NumberRange) } - const observations = yValues.map(getRawValue) as number[] - const numericX = xValues.map(getRawValue) as number[] + const observations = yValues.map(getRawValue) + const numericX = xValues.map(getRawValue) const predictors = observations.map((_, i) => { if (shape.orientation === 'paired') { return [numericX[i]] From a6ba94e5bd666b683594edc5d6d877ba94e30c85 Mon Sep 17 00:00:00 2001 From: tobiadefami Date: Fri, 11 Sep 2026 21:27:08 +0100 Subject: [PATCH 4/6] fix(HF-221): stabilize LINEST standard errors --- docs/guide/known-limitations.md | 2 ++ docs/guide/list-of-differences.md | 2 ++ .../plugin/regression/LinearRegression.ts | 16 +++++++++------- 3 files changed, 13 insertions(+), 7 deletions(-) diff --git a/docs/guide/known-limitations.md b/docs/guide/known-limitations.md index 3a6ceb1e4c..ba47769d71 100644 --- a/docs/guide/known-limitations.md +++ b/docs/guide/known-limitations.md @@ -69,6 +69,8 @@ a circular reference. * If an input expression returns a different predictor count than predicted, `LINEST` returns `#VALUE!`. For example, filtering predictor columns can change the required result width. +* Very large or very small input magnitudes can reduce numerical accuracy. Intermediate calculations, such as residual sums of squares, can still overflow or underflow, producing `#NUM!` or inaccurate statistics. Rescale inputs to more moderate units where possible. See [LINEST numerical differences](list-of-differences.md#linest) for compatibility considerations. + ### OFFSET function HyperFormula resolves the OFFSET function at parse time rather than during evaluation. The parser inspects the arguments and rewrites the expression into a plain cell reference or range. This keeps the dependency graph accurate but imposes several restrictions. diff --git a/docs/guide/list-of-differences.md b/docs/guide/list-of-differences.md index f53af063dc..e4a3b28fdc 100644 --- a/docs/guide/list-of-differences.md +++ b/docs/guide/list-of-differences.md @@ -40,6 +40,8 @@ Unlike Excel, HyperFormula requires a constant `stats` argument because the outp Numerical results can differ for nearly dependent predictors and nearly perfect fits. Coefficient standard errors are sensitive to conditioning, and the F statistic is sensitive to residuals close to machine precision. For an effectively perfect multiple regression, Excel and HyperFormula can return different large finite F values even when the coefficients agree. +Very large or very small input scales can also cause substantial differences in coefficients and statistics, including for a single predictor with a non-perfect fit. HyperFormula does not reproduce Excel's loss of predictors or zero standard errors observed at extreme scales. Rescaling inputs to more moderate units can reduce numerical errors in both engines. + ## General functionalities | Functionality | Examples | HyperFormula | Google Sheets | Microsoft Excel | diff --git a/src/interpreter/plugin/regression/LinearRegression.ts b/src/interpreter/plugin/regression/LinearRegression.ts index d3a286db0b..7b9d793927 100644 --- a/src/interpreter/plugin/regression/LinearRegression.ts +++ b/src/interpreter/plugin/regression/LinearRegression.ts @@ -172,10 +172,12 @@ export function fitLinearRegression(predictors: number[][], observations: number const degreesOfFreedom = n - count const variance = degreesOfFreedom === 0 ? 0 : residualSumSquares / degreesOfFreedom const standardErrors = Array(k).fill(0) - let interceptVariance = 0 + let interceptError = 0 if (statistics && variance !== 0) { - // Each solve produces one column of R^-1. Accumulate only covariance diagonals - // and the intercept's quadratic form, avoiding a full covariance matrix. + const residualStandardError = Math.sqrt(variance) + // Each solve produces one column of R^-1 without a full covariance matrix. + // Accumulate scaled uncertainty magnitudes with hypot: squaring first can + // overflow or underflow even when the final standard error is representable. for (let j = 0; j < count; j++) { const right = Array(count).fill(0) right[j] = 1 @@ -183,17 +185,17 @@ export function fitLinearRegression(predictors: number[][], observations: number let interceptWeight = fitIntercept ? inverseColumn[0] : 0 for (let i = offset; i < count; i++) { const original = columns[i].index - standardErrors[original] += variance * inverseColumn[i] ** 2 + standardErrors[original] = Math.hypot(standardErrors[original], residualStandardError * inverseColumn[i]) interceptWeight -= means[original] * inverseColumn[i] } - interceptVariance += variance * interceptWeight ** 2 + interceptError = Math.hypot(interceptError, residualStandardError * interceptWeight) } } return { coefficients, intercept, - standardErrors: standardErrors.map(Math.sqrt), - interceptError: Math.sqrt(interceptVariance), + standardErrors, + interceptError, residualSumSquares, totalSumSquares, degreesOfFreedom, From b8a20aa0aff87a38cc4e53b1fe28656e57f5d96f Mon Sep 17 00:00:00 2001 From: tobiadefami Date: Mon, 14 Sep 2026 19:57:01 +0100 Subject: [PATCH 5/6] refactor(HF-221): clarify LINEST regression flow --- src/interpreter/plugin/RegressionPlugin.ts | 45 ++- .../plugin/regression/LinearRegression.ts | 314 ++++++++++++------ 2 files changed, 241 insertions(+), 118 deletions(-) diff --git a/src/interpreter/plugin/RegressionPlugin.ts b/src/interpreter/plugin/RegressionPlugin.ts index 7810eb2e7c..9df3160c1e 100644 --- a/src/interpreter/plugin/RegressionPlugin.ts +++ b/src/interpreter/plugin/RegressionPlugin.ts @@ -16,24 +16,46 @@ import {fitLinearRegression, LinearRegressionResult} from './regression/LinearRe /** The orientation and predictor count shared by size prediction and runtime validation. */ interface RegressionShape { - predictors: number, + predictorCount: number, orientation: 'paired' | 'columns' | 'rows', } /** Classifies the original dimensions before any values are flattened. */ function regressionShape(y: ArraySize, x?: ArraySize): RegressionShape | undefined { if (x === undefined || (y.width === x.width && y.height === x.height)) { - return {predictors: 1, orientation: 'paired'} + return {predictorCount: 1, orientation: 'paired'} } if (y.width === 1 && y.height === x.height) { - return {predictors: x.width, orientation: 'columns'} + return {predictorCount: x.width, orientation: 'columns'} } if (y.height === 1 && y.width === x.width) { - return {predictors: x.height, orientation: 'rows'} + return {predictorCount: x.height, orientation: 'rows'} } return undefined } +/** + * Converts validated, flattened predictor values into one row per observation. + * In columns orientation, each source row is an observation; in rows orientation, + * each source row is a predictor. Paired ranges supply one predictor per observation. + */ +function buildPredictorRows(predictorValues: number[], observationCount: number, shape: RegressionShape): number[][] { + return Array.from({length: observationCount}, (_, observationIndex) => { + if (shape.orientation === 'paired') { + return [predictorValues[observationIndex]] + } + return Array.from({length: shape.predictorCount}, (_, predictorIndex) => { + let flatIndex: number + if (shape.orientation === 'columns') { + flatIndex = observationIndex * shape.predictorCount + predictorIndex + } else { + flatIndex = predictorIndex * observationCount + observationIndex + } + return predictorValues[flatIndex] + }) + }) +} + /** LINEST rejects empty strings as Boolean options, including formula-generated strings. */ function regressionBoolean(value: InternalScalarValue): boolean | CellError | undefined { return value === '' ? undefined : coerceScalarToBoolean(value) @@ -154,10 +176,10 @@ export class RegressionPlugin extends FunctionPlugin implements FunctionPluginTy if (shape === undefined) { return new CellError(ErrorType.REF, ErrorMessage.ArrayDimensions) } - if (shape.predictors + 1 > this.config.maxColumns) { + if (shape.predictorCount + 1 > this.config.maxColumns) { return new CellError(ErrorType.VALUE, ErrorMessage.ValueLarge) } - if (this.linestArraySize(ast, state).width !== shape.predictors + 1) { + if (this.linestArraySize(ast, state).width !== shape.predictorCount + 1) { return new CellError(ErrorType.VALUE, ErrorMessage.LinestStaticSize) } const yValues = Array.from(knownY.valuesFromTopLeftCorner()) @@ -167,12 +189,7 @@ export class RegressionPlugin extends FunctionPlugin implements FunctionPluginTy } const observations = yValues.map(getRawValue) const numericX = xValues.map(getRawValue) - const predictors = observations.map((_, i) => { - if (shape.orientation === 'paired') { - return [numericX[i]] - } - return Array.from({length: shape.predictors}, (_, j) => numericX[shape.orientation === 'columns' ? i * shape.predictors + j : j * observations.length + i]) - }) + const predictors = buildPredictorRows(numericX, observations.length, shape) const fit = fitLinearRegression(predictors, observations, fitIntercept, statistics) return regressionOutput(fit, statistics, fitIntercept) }) @@ -191,9 +208,9 @@ export class RegressionPlugin extends FunctionPlugin implements FunctionPluginTy const y = this.arraySizeForAst(ast.args[0], arrayState) const x = ast.args.length < 2 || ast.args[1].type === AstNodeType.EMPTY ? undefined : this.arraySizeForAst(ast.args[1], arrayState) const shape = regressionShape(y, x) - if (shape === undefined || shape.predictors + 1 > this.config.maxColumns) { + if (shape === undefined || shape.predictorCount + 1 > this.config.maxColumns) { return ArraySize.error() } - return new ArraySize(shape.predictors + 1, statistics ? 5 : 1) + return new ArraySize(shape.predictorCount + 1, statistics ? 5 : 1) } } diff --git a/src/interpreter/plugin/regression/LinearRegression.ts b/src/interpreter/plugin/regression/LinearRegression.ts index 7b9d793927..8cdf03b019 100644 --- a/src/interpreter/plugin/regression/LinearRegression.ts +++ b/src/interpreter/plugin/regression/LinearRegression.ts @@ -15,13 +15,31 @@ export interface LinearRegressionResult { retainedPredictorCount: number, } -/** A predictor column with its original position and pre-factorization norm. */ +/** + * Column-major storage: columns[columnIndex].values[rowIndex] is a matrix entry. + * Factorization overwrites the values with the triangular system. predictorIndex + * preserves the original predictor position through swaps; -1 identifies the intercept. + */ interface RegressionColumn { values: number[], - index: number, + predictorIndex: number, + /** Norm after optional centering, before any reflections. */ originalNorm: number, } +/** A reflection and the leading value it produces in the selected column tail. */ +interface HouseholderReflection { + vector: number[], + factor: number, + diagonalValue: number, +} + +/** Coefficient uncertainties restored to the original predictor order. */ +interface CoefficientStandardErrors { + standardErrors: number[], + interceptError: number, +} + /** * Rank threshold relative to a predictor's original centered norm. * Excel Online retains a 1e-5 perturbation of a duplicated predictor and removes 1e-6. @@ -30,10 +48,10 @@ interface RegressionColumn { const RANK_TOLERANCE = 1e-6 /** Computes a Euclidean norm without squaring large inputs or spreading an unbounded range. */ -function norm(values: number[], start = 0): number { +function norm(values: number[], startRow = 0): number { let result = 0 - for (let i = start; i < values.length; i++) { - result = Math.hypot(result, values[i]) + for (let rowIndex = startRow; rowIndex < values.length; rowIndex++) { + result = Math.hypot(result, values[rowIndex]) } return result } @@ -48,42 +66,143 @@ function mean(values: number[]): number { return origin + sum } -/** Applies a Householder reflection to the trailing part of a column. */ -function reflect(values: number[], vector: number[], start: number, factor: number): void { - let product = 0 - for (let i = start; i < values.length; i++) { - product += vector[i - start] * values[i] +/** + * Constructs a reflection for a nonzero column tail without modifying the column. + * + * Let x = values.slice(startRow), r = magnitude, and e0 = [1, 0, ...]. + * Choose s = 1 when x[0] >= 0, otherwise -1. + * The vector v = x / r + s * e0 reflects x to [-s * r, 0, ...]. Adding the same + * sign avoids cancellation when constructing v's first entry. + * Its squared length is v^T * v = 2 * (1 + abs(x[0]) / r), so the reflection + * factor 2 / (v^T * v) simplifies to 1 / (1 + abs(x[0]) / r). + * + * @param values - selected column, including any already processed rows + * @param startRow - first row of the column tail to reflect + * @param magnitude - nonzero norm of values starting at startRow, already computed by the caller + */ +function createHouseholderReflection(values: number[], startRow: number, magnitude: number): HouseholderReflection { + const sign = values[startRow] >= 0 ? 1 : -1 + const vector = values.slice(startRow).map(value => value / magnitude) + vector[0] += sign + const factor = 1 / (1 + Math.abs(values[startRow]) / magnitude) + return {vector, factor, diagonalValue: -sign * magnitude} +} + +/** + * Applies a Householder reflection to values starting at startRow, in place. + * + * For the active tail w and reflection vector v, with ^T denoting transpose: + * dotProduct = v^T * w + * reflectionScale = factor * dotProduct, where factor = 2 / (v^T * v) + * w_new = w - reflectionScale * v + * The subtracted vector is twice the projection onto v: this reverses that + * component while preserving the perpendicular component and the total length. + * + * startRow stays fixed during this call; rowIndex - startRow indexes the shorter v. + * Earlier rows remain unchanged. Compute the full dot product before overwriting w. + */ +function applyHouseholderReflection(values: number[], vector: number[], startRow: number, factor: number): void { + let dotProduct = 0 + for (let rowIndex = startRow; rowIndex < values.length; rowIndex++) { + dotProduct += vector[rowIndex - startRow] * values[rowIndex] + } + const reflectionScale = dotProduct * factor + for (let rowIndex = startRow; rowIndex < values.length; rowIndex++) { + values[rowIndex] -= reflectionScale * vector[rowIndex - startRow] } - product *= factor - for (let i = start; i < values.length; i++) { - values[i] -= product * vector[i - start] +} + +/** + * Selects the largest remaining column tail without changing column order. + * diagonalIndex is both the first candidate column and the first active row. + * Strict comparison keeps the first candidate on ties, preserving Excel's duplicate selection. + */ +function selectPivotColumn(columns: RegressionColumn[], diagonalIndex: number, activeColumnCount: number): number { + let selectedColumnIndex = diagonalIndex + let largestNorm = -1 + for (let columnIndex = diagonalIndex; columnIndex < activeColumnCount; columnIndex++) { + const trailingNorm = norm(columns[columnIndex].values, diagonalIndex) + if (trailingNorm > largestNorm) { + largestNorm = trailingNorm + selectedColumnIndex = columnIndex + } } + return selectedColumnIndex } -/** Solves the retained triangular system, treating an exactly zero predictor as a zero coefficient. */ -function backSubstitute(columns: RegressionColumn[], count: number, right: number[]): number[] { - const solution = Array(count).fill(0) - for (let i = count - 1; i >= 0; i--) { - if (columns[i].values[i] === 0) { +/** + * Solves R * solution = rightHandSide from the bottom row upward. + * R[row, column] is stored as columns[column].values[row]. Subtract the already + * solved terms to the right of the diagonal, then divide by the diagonal value. + * Exactly zero diagonals leave a zero coefficient for Excel-compatible zero columns. + */ +function backSubstitute(columns: RegressionColumn[], processedColumnCount: number, rightHandSide: number[]): number[] { + const solution = Array(processedColumnCount).fill(0) + for (let rowIndex = processedColumnCount - 1; rowIndex >= 0; rowIndex--) { + if (columns[rowIndex].values[rowIndex] === 0) { continue } - let value = right[i] - for (let j = i + 1; j < count; j++) { - value -= columns[j].values[i] * solution[j] + let remainingValue = rightHandSide[rowIndex] + for (let columnIndex = rowIndex + 1; columnIndex < processedColumnCount; columnIndex++) { + remainingValue -= columns[columnIndex].values[rowIndex] * solution[columnIndex] } - solution[i] = value / columns[i].values[i] + solution[rowIndex] = remainingValue / columns[rowIndex].values[rowIndex] } return solution } +/** + * Computes standard errors from the retained triangular system, in original predictor order. + * + * For nonsingular R, covariance = residualVariance * inverseR * transpose(inverseR). + * A coefficient's standard error is therefore sqrt(residualVariance) times the norm + * of its row of inverseR. Solve R * inverseColumn = e_j for each unit vector e_j + * to accumulate those row norms without building a full inverse or covariance matrix. + * + * With an intercept, centering changes it by -sum(coefficient * predictorMean), so each + * inverse column contributes interceptWeight = inverseColumn[0] - sum(mean * inverseEntry). + * Zero diagonals retain backSubstitute's Excel-compatible zero-coefficient behavior. + */ +function computeCoefficientStandardErrors( + columns: RegressionColumn[], + processedColumnCount: number, + predictorMeans: number[], + fitIntercept: boolean, + residualVariance: number, +): CoefficientStandardErrors { + const standardErrors = Array(predictorMeans.length).fill(0) + let interceptError = 0 + if (residualVariance === 0) { + return {standardErrors, interceptError} + } + + const interceptColumnCount = fitIntercept ? 1 : 0 + const residualStandardError = Math.sqrt(residualVariance) + for (let inverseColumnIndex = 0; inverseColumnIndex < processedColumnCount; inverseColumnIndex++) { + const rightHandSide = Array(processedColumnCount).fill(0) + rightHandSide[inverseColumnIndex] = 1 + const inverseColumn = backSubstitute(columns, processedColumnCount, rightHandSide) + let interceptWeight = fitIntercept ? inverseColumn[0] : 0 + for (let columnIndex = interceptColumnCount; columnIndex < processedColumnCount; columnIndex++) { + const predictorIndex = columns[columnIndex].predictorIndex + // Scale before hypot: squaring inverse entries first can overflow or underflow. + standardErrors[predictorIndex] = Math.hypot(standardErrors[predictorIndex], residualStandardError * inverseColumn[columnIndex]) + interceptWeight -= predictorMeans[predictorIndex] * inverseColumn[columnIndex] + } + interceptError = Math.hypot(interceptError, residualStandardError * interceptWeight) + } + return {standardErrors, interceptError} +} + /** * Fits observations using pivoted Householder QR, with a fixed leading intercept column. - * Predictor centering preserves accuracy for large offsets. The factorization is column-major, - * uses O(n*k + k*k) storage, and costs O(n*k*k) for n >= k; Q is never materialized. + * Prepares columns, factorizes them, restores coefficients, then computes fit statistics. + * Predictor centering preserves accuracy for large offsets. For n observations and k predictors, + * the factorization uses O(n*k + k*k) storage and costs O(n*k*k) for n >= k; Q is never materialized. * * Exactly zero predictors intentionally consume an available QR slot. Excel Online returns a * zero coefficient but still excludes that transformed response entry from residual statistics. - * A dependent nonzero predictor is instead removed using its relative norm, increasing df. + * A dependent nonzero predictor is instead removed using its relative norm, increasing degrees of freedom. * * @param predictors - observation rows, with predictor columns in their original order * @param observations - numeric response for each observation @@ -91,106 +210,93 @@ function backSubstitute(columns: RegressionColumn[], count: number, right: numbe * @param statistics - whether to compute coefficient uncertainties */ export function fitLinearRegression(predictors: number[][], observations: number[], fitIntercept: boolean, statistics: boolean): LinearRegressionResult { - const n = observations.length - const k = predictors[0].length - const offset = fitIntercept ? 1 : 0 - const means = Array(k).fill(0) + const observationCount = observations.length + const predictorCount = predictors[0].length + const interceptColumnCount = fitIntercept ? 1 : 0 + const predictorMeans = Array(predictorCount).fill(0) const columns: RegressionColumn[] = [] if (fitIntercept) { - columns.push({values: Array(n).fill(1), index: -1, originalNorm: Math.sqrt(n)}) + columns.push({values: Array(observationCount).fill(1), predictorIndex: -1, originalNorm: Math.sqrt(observationCount)}) } - for (let j = 0; j < k; j++) { - const values = predictors.map(row => row[j]) - means[j] = fitIntercept ? mean(values) : 0 - const centered = values.map(value => value - means[j]) - columns.push({values: centered, index: j, originalNorm: norm(centered)}) + for (let predictorIndex = 0; predictorIndex < predictorCount; predictorIndex++) { + const values = predictors.map(row => row[predictorIndex]) + predictorMeans[predictorIndex] = fitIntercept ? mean(values) : 0 + const centeredValues = values.map(value => value - predictorMeans[predictorIndex]) + columns.push({values: centeredValues, predictorIndex, originalNorm: norm(centeredValues)}) } - // Treat the intercept as the last input column moved to the front by a swap. - // Retain this order when pivot norms tie, so duplicated predictors agree with Excel. - if (fitIntercept && k > 1) { + // Excel tie ordering: [intercept, x1, x2, ...] -> [intercept, x2, ..., x1]. + // This matches swapping a trailing intercept into the first input position. + if (fitIntercept && predictorCount > 1) { const firstPredictor = columns.splice(1, 1)[0] columns.push(firstPredictor) } // Remove a common response offset before reflection to preserve small variations. const responseOrigin = fitIntercept ? observations[0] : 0 - const transformedY = observations.map(value => value - responseOrigin) - let remaining = columns.length - let count = 0 - while (count < Math.min(n, remaining)) { - if (count >= offset) { - let selected = count - let largest = -1 - for (let j = count; j < remaining; j++) { - const magnitude = norm(columns[j].values, count) - if (magnitude > largest) { - largest = magnitude - selected = j - } - } - [columns[count], columns[selected]] = [columns[selected], columns[count]] + const transformedResponse = observations.map(value => value - responseOrigin) + + // Columns are partitioned as [processed | candidates | rejected]: + // [0, processedColumnCount), [processedColumnCount, activeColumnCount), and the rest. + // processedColumnCount is also the current diagonal index, not the mathematical rank: + // exactly zero columns still consume a slot for Excel-compatible residual statistics. + let activeColumnCount = columns.length + let processedColumnCount = 0 + while (processedColumnCount < Math.min(observationCount, activeColumnCount)) { + if (processedColumnCount >= interceptColumnCount) { + const pivotColumnIndex = selectPivotColumn(columns, processedColumnCount, activeColumnCount) + const previousColumn = columns[processedColumnCount] + columns[processedColumnCount] = columns[pivotColumnIndex] + columns[pivotColumnIndex] = previousColumn } - const column = columns[count] - const magnitude = norm(column.values, count) - if (count >= offset && magnitude < RANK_TOLERANCE * column.originalNorm) { - remaining-- - columns[count] = columns[remaining] - columns[remaining] = column + const column = columns[processedColumnCount] + const trailingNorm = norm(column.values, processedColumnCount) + // Rejected predictors do not advance the diagonal. An originally zero column + // gives 0 < 0 here, so it is retained and advances processedColumnCount below. + if (processedColumnCount >= interceptColumnCount && trailingNorm < RANK_TOLERANCE * column.originalNorm) { + activeColumnCount-- + columns[processedColumnCount] = columns[activeColumnCount] + columns[activeColumnCount] = column continue } - if (magnitude !== 0) { - const sign = column.values[count] >= 0 ? 1 : -1 - const vector = column.values.slice(count).map(value => value / magnitude) - vector[0] += sign - const factor = 1 / (1 + Math.abs(column.values[count]) / magnitude) - for (let j = count + 1; j < remaining; j++) { - reflect(columns[j].values, vector, count, factor) + if (trailingNorm !== 0) { + const {vector, factor, diagonalValue} = createHouseholderReflection(column.values, processedColumnCount, trailingNorm) + for (let columnIndex = processedColumnCount + 1; columnIndex < activeColumnCount; columnIndex++) { + applyHouseholderReflection(columns[columnIndex].values, vector, processedColumnCount, factor) } - reflect(transformedY, vector, count, factor) - column.values[count] = -sign * magnitude - column.values.fill(0, count + 1) + applyHouseholderReflection(transformedResponse, vector, processedColumnCount, factor) + column.values[processedColumnCount] = diagonalValue + column.values.fill(0, processedColumnCount + 1) } - count++ + processedColumnCount++ } - const solution = backSubstitute(columns, count, transformedY) - const coefficients = Array(k).fill(0) - for (let j = offset; j < count; j++) { - coefficients[columns[j].index] = solution[j] + // Solve in pivot order, then restore the caller's predictor order. Rejected slots stay zero. + const solution = backSubstitute(columns, processedColumnCount, transformedResponse) + const coefficients = Array(predictorCount).fill(0) + for (let columnIndex = interceptColumnCount; columnIndex < processedColumnCount; columnIndex++) { + coefficients[columns[columnIndex].predictorIndex] = solution[columnIndex] } - const intercept = fitIntercept ? mean(observations) - coefficients.reduce((sum, value, j) => sum + value * means[j], 0) : 0 - const yMean = fitIntercept ? mean(observations) : 0 - const centeredY = observations.map(value => value - yMean) - const totalSumSquares = norm(centeredY) ** 2 - let residualSumSquares = norm(transformedY, count) ** 2 + // Recover the intercept in original coordinates: b = mean(y) - sum(coefficient * mean(x)). + const intercept = fitIntercept + ? mean(observations) - coefficients.reduce((sum, coefficient, predictorIndex) => sum + coefficient * predictorMeans[predictorIndex], 0) + : 0 + + const responseMean = fitIntercept ? mean(observations) : 0 + const centeredResponse = observations.map(value => value - responseMean) + const totalSumSquares = norm(centeredResponse) ** 2 + // Reflections preserve squared error; entries below the solved rows are the residual tail. + let residualSumSquares = norm(transformedResponse, processedColumnCount) ** 2 // Simple exact fits have zero residual statistics in Excel. Discard only // reflection roundoff, whose squared scale is O((n * machine epsilon)^2). - const simpleFitRoundoff = k === 1 && residualSumSquares <= totalSumSquares * (n * Number.EPSILON) ** 2 - if (totalSumSquares === 0 || count === n || simpleFitRoundoff) { + const simpleFitRoundoff = predictorCount === 1 && residualSumSquares <= totalSumSquares * (observationCount * Number.EPSILON) ** 2 + if (totalSumSquares === 0 || processedColumnCount === observationCount || simpleFitRoundoff) { residualSumSquares = 0 } - const degreesOfFreedom = n - count - const variance = degreesOfFreedom === 0 ? 0 : residualSumSquares / degreesOfFreedom - const standardErrors = Array(k).fill(0) - let interceptError = 0 - if (statistics && variance !== 0) { - const residualStandardError = Math.sqrt(variance) - // Each solve produces one column of R^-1 without a full covariance matrix. - // Accumulate scaled uncertainty magnitudes with hypot: squaring first can - // overflow or underflow even when the final standard error is representable. - for (let j = 0; j < count; j++) { - const right = Array(count).fill(0) - right[j] = 1 - const inverseColumn = backSubstitute(columns, count, right) - let interceptWeight = fitIntercept ? inverseColumn[0] : 0 - for (let i = offset; i < count; i++) { - const original = columns[i].index - standardErrors[original] = Math.hypot(standardErrors[original], residualStandardError * inverseColumn[i]) - interceptWeight -= means[original] * inverseColumn[i] - } - interceptError = Math.hypot(interceptError, residualStandardError * interceptWeight) - } - } + const degreesOfFreedom = observationCount - processedColumnCount + const residualVariance = degreesOfFreedom === 0 ? 0 : residualSumSquares / degreesOfFreedom + const {standardErrors, interceptError} = statistics + ? computeCoefficientStandardErrors(columns, processedColumnCount, predictorMeans, fitIntercept, residualVariance) + : {standardErrors: Array(predictorCount).fill(0), interceptError: 0} return { coefficients, intercept, @@ -199,6 +305,6 @@ export function fitLinearRegression(predictors: number[][], observations: number residualSumSquares, totalSumSquares, degreesOfFreedom, - retainedPredictorCount: count - offset, + retainedPredictorCount: processedColumnCount - interceptColumnCount, } } From 558bd2892edab46eeffe2243dccc72397fbf463a Mon Sep 17 00:00:00 2001 From: tobiadefami Date: Wed, 16 Sep 2026 21:07:11 +0100 Subject: [PATCH 6/6] refactor(HF-221): separate LINEST calculation stages --- .../plugin/regression/LinearRegression.ts | 126 ++++++++++++++---- 1 file changed, 101 insertions(+), 25 deletions(-) diff --git a/src/interpreter/plugin/regression/LinearRegression.ts b/src/interpreter/plugin/regression/LinearRegression.ts index 8cdf03b019..86b7b1400a 100644 --- a/src/interpreter/plugin/regression/LinearRegression.ts +++ b/src/interpreter/plugin/regression/LinearRegression.ts @@ -27,6 +27,21 @@ interface RegressionColumn { originalNorm: number, } +/** Working arrays owned by one fit, with centering offsets kept in the original predictor order. */ +interface PreparedRegression { + columns: RegressionColumn[], + predictorMeans: number[], + transformedResponse: number[], +} + +/** Residual statistics shared by coefficient uncertainty calculation and spreadsheet output. */ +interface ResidualStatistics { + totalSumSquares: number, + residualSumSquares: number, + degreesOfFreedom: number, + residualVariance: number, +} + /** A reflection and the leading value it produces in the selected column tail. */ interface HouseholderReflection { vector: number[], @@ -42,7 +57,7 @@ interface CoefficientStandardErrors { /** * Rank threshold relative to a predictor's original centered norm. - * Excel Online retains a 1e-5 perturbation of a duplicated predictor and removes 1e-6. + * Sampled Excel Online cases with proportional predictors retain a 1e-5 perturbation and remove 1e-6. * Keep the threshold separate from roundoff handling for the fitted statistics. */ const RANK_TOLERANCE = 1e-6 @@ -195,24 +210,20 @@ function computeCoefficientStandardErrors( } /** - * Fits observations using pivoted Householder QR, with a fixed leading intercept column. - * Prepares columns, factorizes them, restores coefficients, then computes fit statistics. - * Predictor centering preserves accuracy for large offsets. For n observations and k predictors, - * the factorization uses O(n*k + k*k) storage and costs O(n*k*k) for n >= k; Q is never materialized. + * Copies observations and predictors into working arrays without modifying the inputs. + * With an intercept, uses centeredX = x - mean(x) and shiftedY = y - y[0] to preserve + * small variations on large offsets. Without an intercept, leaves both coordinates unchanged. + * Predictor means retain their input order for recovering b = mean(y) - sum(m * mean(x)). + * When no intercept is fitted, predictorMeans contains zero offsets instead. * - * Exactly zero predictors intentionally consume an available QR slot. Excel Online returns a - * zero coefficient but still excludes that transformed response entry from residual statistics. - * A dependent nonzero predictor is instead removed using its relative norm, increasing degrees of freedom. - * - * @param predictors - observation rows, with predictor columns in their original order + * @param predictors - observation rows with predictor columns in their original order * @param observations - numeric response for each observation - * @param fitIntercept - whether to include a constant term - * @param statistics - whether to compute coefficient uncertainties + * @param fitIntercept - whether to add a column of ones and shift the coordinates + * @returns Working columns in Excel tie order, predictor centering offsets, and the shifted response. */ -export function fitLinearRegression(predictors: number[][], observations: number[], fitIntercept: boolean, statistics: boolean): LinearRegressionResult { +function prepareRegression(predictors: number[][], observations: number[], fitIntercept: boolean): PreparedRegression { const observationCount = observations.length const predictorCount = predictors[0].length - const interceptColumnCount = fitIntercept ? 1 : 0 const predictorMeans = Array(predictorCount).fill(0) const columns: RegressionColumn[] = [] if (fitIntercept) { @@ -234,7 +245,26 @@ export function fitLinearRegression(predictors: number[][], observations: number // Remove a common response offset before reflection to preserve small variations. const responseOrigin = fitIntercept ? observations[0] : 0 const transformedResponse = observations.map(value => value - responseOrigin) + return {columns, predictorMeans, transformedResponse} +} +/** + * Overwrites working columns with a pivoted triangular system and applies the same + * Householder reflections to the response. Retained columns store R; the transformed + * response supplies the right-hand side of R * solution = response and its residual tail. + * Original predictor positions remain available through each column's predictorIndex. + * + * Exactly zero predictors intentionally consume an available slot. Excel Online returns + * a zero coefficient but still excludes that response entry from residual statistics. + * A dependent nonzero predictor is removed without consuming a slot. + * + * @param columns - mutable working columns, including a fixed leading intercept when present + * @param transformedResponse - mutable shifted response, transformed in place alongside the columns + * @param interceptColumnCount - one when fitting an intercept, otherwise zero + * @returns Number of processed columns, including zero columns; this is not mathematical rank. + */ +function factorizeRegressionInPlace(columns: RegressionColumn[], transformedResponse: number[], interceptColumnCount: number): number { + const observationCount = transformedResponse.length // Columns are partitioned as [processed | candidates | rejected]: // [0, processedColumnCount), [processedColumnCount, activeColumnCount), and the rest. // processedColumnCount is also the current diagonal index, not the mathematical rank: @@ -269,18 +299,30 @@ export function fitLinearRegression(predictors: number[][], observations: number } processedColumnCount++ } + return processedColumnCount +} - // Solve in pivot order, then restore the caller's predictor order. Rejected slots stay zero. - const solution = backSubstitute(columns, processedColumnCount, transformedResponse) - const coefficients = Array(predictorCount).fill(0) - for (let columnIndex = interceptColumnCount; columnIndex < processedColumnCount; columnIndex++) { - coefficients[columns[columnIndex].predictorIndex] = solution[columnIndex] - } - // Recover the intercept in original coordinates: b = mean(y) - sum(coefficient * mean(x)). - const intercept = fitIntercept - ? mean(observations) - coefficients.reduce((sum, coefficient, predictorIndex) => sum + coefficient * predictorMeans[predictorIndex], 0) - : 0 - +/** + * Computes response variation and residual error without modifying either response array. + * Total variation is sum((y - mean(y))^2) with an intercept, otherwise sum(y^2). + * Residual error is the squared norm below the processed rows of the transformed response. + * Uses the factorization's Excel-compatible row boundary, including consumed zero columns. + * + * @param observations - response values in their original coordinates + * @param transformedResponse - shifted response after the factorization's reflections + * @param processedColumnCount - number of response entries assigned to the triangular solve + * @param predictorCount - original predictor count, used for the simple-fit roundoff rule + * @param fitIntercept - whether total variation is measured around the response mean + * @returns Sums of squares, remaining degrees of freedom, and residual variance (zero when no degrees remain). + */ +function computeResidualStatistics( + observations: number[], + transformedResponse: number[], + processedColumnCount: number, + predictorCount: number, + fitIntercept: boolean, +): ResidualStatistics { + const observationCount = observations.length const responseMean = fitIntercept ? mean(observations) : 0 const centeredResponse = observations.map(value => value - responseMean) const totalSumSquares = norm(centeredResponse) ** 2 @@ -294,6 +336,40 @@ export function fitLinearRegression(predictors: number[][], observations: number } const degreesOfFreedom = observationCount - processedColumnCount const residualVariance = degreesOfFreedom === 0 ? 0 : residualSumSquares / degreesOfFreedom + return {totalSumSquares, residualSumSquares, degreesOfFreedom, residualVariance} +} + +/** + * Fits observations using pivoted Householder QR, with a fixed leading intercept column. + * Prepares columns, factorizes them, restores coefficients, then computes fit statistics. + * Predictor centering preserves accuracy for large offsets. For n observations and k predictors, + * the factorization uses O(n*k + k*k) storage and costs O(n*k*k) for n >= k; Q is never materialized. + * + * @param predictors - observation rows, with predictor columns in their original order + * @param observations - numeric response for each observation + * @param fitIntercept - whether to include a constant term + * @param statistics - whether to compute coefficient uncertainties + */ +export function fitLinearRegression(predictors: number[][], observations: number[], fitIntercept: boolean, statistics: boolean): LinearRegressionResult { + const predictorCount = predictors[0].length + const interceptColumnCount = fitIntercept ? 1 : 0 + const {columns, predictorMeans, transformedResponse} = prepareRegression(predictors, observations, fitIntercept) + const processedColumnCount = factorizeRegressionInPlace(columns, transformedResponse, interceptColumnCount) + + // Solve in pivot order, then restore the caller's predictor order. Rejected slots stay zero. + const solution = backSubstitute(columns, processedColumnCount, transformedResponse) + const coefficients = Array(predictorCount).fill(0) + for (let columnIndex = interceptColumnCount; columnIndex < processedColumnCount; columnIndex++) { + coefficients[columns[columnIndex].predictorIndex] = solution[columnIndex] + } + // Recover the intercept in original coordinates: b = mean(y) - sum(coefficient * mean(x)). + const intercept = fitIntercept + ? mean(observations) - coefficients.reduce((sum, coefficient, predictorIndex) => sum + coefficient * predictorMeans[predictorIndex], 0) + : 0 + + const {totalSumSquares, residualSumSquares, degreesOfFreedom, residualVariance} = computeResidualStatistics( + observations, transformedResponse, processedColumnCount, predictorCount, fitIntercept, + ) const {standardErrors, interceptError} = statistics ? computeCoefficientStandardErrors(columns, processedColumnCount, predictorMeans, fitIntercept, residualVariance) : {standardErrors: Array(predictorCount).fill(0), interceptError: 0}