@geodetic-lse/core

Generic linear least squares solver with internal reliability analysis, based on the method of parametric adjustment (indirect observations).

Installation

npm install @geodetic-lse/core

Overview

This package provides the mathematical core of the Geodetic LSE engine:

• Normal equations solver — using LU decomposition for numerical stability
• Reliability analysis — Baarda's normalized residuals, local redundancy numbers, and estimated gross errors

It is designed to be discipline-agnostic: any problem that can be expressed as A·dx + v = dl can use this solver.

API

solveLeastSquares(A, P, l, fx, x)

Solves the weighted least squares system using the normal equations method.

| Parameter | Type | Description | | --------- | -------- | ------------------------------------------------ | | A | Matrix | Design matrix (m × n) | | P | Matrix | Weight matrix (m × m), diagonal | | l | Matrix | Observation vector (m × 1) | | fx | Matrix | Computed values at initial approximation (m × 1) | | x | Matrix | Initial approximation of unknowns (n × 1) |

Returns a LeastSquaresResult:

| Field | Type | Description | | ------------ | -------- | ------------------------------- | | x | Matrix | Adjusted parameters (n × 1) | | v | Matrix | Residuals (m × 1) | | N | Matrix | Normal matrix AᵀPA (n × n) | | Qxx | Matrix | Cofactor matrix N⁻¹ (n × n) | | s0 | number | A posteriori standard deviation | | redundancy | number | Degrees of freedom (m - n) |

computeReliability(sigma0, Qll, P, A, result, alpha?, beta?)

Computes internal and external reliability indicators following Baarda's theory.

| Parameter | Type | Default | Description | | --------- | -------------------- | ------- | --------------------------------------- | | sigma0 | number | | A priori standard deviation | | Qll | Matrix | | Cofactor matrix of observations (m × m) | | P | Matrix | | Weight matrix (m × m) | | A | Matrix | | Design matrix (m × n) | | result | LeastSquaresResult | | Output of solveLeastSquares | | alpha | number | 0.01 | Significance level (type I error rate) | | beta | number | 0.05 | Type II error rate |

Returns a ReliabilityResult:

| Field | Type | Description | | -------- | -------- | ------------------------------------------------------------------------------------------------------------------- | | Qvv | Matrix | Cofactor matrix of residuals (m × m) | | w | Matrix | Normalized residuals — Baarda's wᵢ = vᵢ / (σ₀·√qvvᵢᵢ) (m × 1) | | z | Matrix | Local redundancy numbers rᵢ = qvvᵢᵢ / qllᵢᵢ, ∑rᵢ = redundancy (m × 1) | | g | Matrix | Estimated gross errors gᵢ = −vᵢ / rᵢ; NaN when rᵢ ≈ 0 (m × 1) | | nablaL | Matrix | Internal reliability ∇lᵢ = δ₀·σ₀·√qllᵢᵢ / √rᵢ — min detectable bias; NaN when rᵢ ≈ 0 (m × 1) | | nablaX | Matrix | External reliability ∇xᵢ = Qxx·Aᵀ·P·eᵢ·∇lᵢ — effect on unknowns per observation; NaN column when ∇lᵢ is NaN (n × m) |

solveAndAnalyze(A, P, l, fx, x, sigma0, Qll)

Convenience wrapper that calls solveLeastSquares then computeReliability with default alpha = 0.01 and beta = 0.05, and returns a merged result object containing all fields from both LeastSquaresResult and ReliabilityResult.

Usage

import { matrix } from 'mathjs'
import { solveLeastSquares, computeReliability } from '@geodetic-lse/core'

// Design matrix, weight matrix, observations, computed values, initial unknowns
const A = matrix([
    [1, 0],
    [0, 1],
    [1, -1],
])
const P = matrix([
    [1, 0, 0],
    [0, 1, 0],
    [0, 0, 1],
])
const l = matrix([[10.01], [20.03], [-10.02]])
const fx = matrix([[10.0], [20.0], [-10.0]])
const x0 = matrix([[10.0], [20.0]])

const result = solveLeastSquares(A, P, l, fx, x0)

console.log('Adjusted parameters:', result.x.toArray())
console.log('Residuals:', result.v.toArray())
console.log('s0:', result.s0)

Mathematical model

The functional model is:

A·dx + v = dl     where  dl = l − f(x₀)

The normal equations are formed and solved via LU decomposition (Doolittle / lup):

N·dx = AᵀP·dl     where  N = AᵀPA
Qxx  = N⁻¹
s0   = sqrt(vᵀPv / r)

An ill-conditioning guard is applied: if the 1-norm condition number κ₁(N) exceeds 10¹², an error is thrown.

Reliability

Internal and external reliability follow Baarda's theory. The non-centrality parameter δ₀ is derived from the chosen significance level α (type I error) and power 1−β (type II error):

δ₀ = Φ⁻¹(1 − α/2) − Φ⁻¹(β)

where Φ⁻¹ is the inverse standard normal CDF. With defaults α = 0.01 and β = 0.05:

Internal reliability:  ∇lᵢ = δ₀·σ₀·√qllᵢᵢ / √rᵢ   (min detectable bias on obs i)
External reliability:  ∇xᵢ = Qxx·Aᵀ·P·eᵢ·∇lᵢ       (effect on unknowns, n×1 per obs i)

∇lᵢ and the i-th column of ∇xᵢ are NaN when rᵢ ≈ 0 (observation has no redundancy).

License

MPL-2.0