Least-Squares Graph Methods

• Least-squares graph is a family of formulations that leverages graph topology to structure the minimization of squared residuals across diverse applications.
• These methods exploit graph Laplacians, factor graphs, and consensus dynamics to achieve enhanced numerical efficiency, robustness, and scalability.
• They integrate classical ranking, sensor network calibration, and online signal estimation by enforcing local constraints and spectral properties within the graph.

Taken together, these works suggest that least-squares graph is best understood as a family of formulations in which least-squares estimation is organized by graph structure rather than as a single algorithm. In different settings, the graph may be a weighted, oriented graph carrying pairwise comparisons on edges, a comparison multigraph whose Laplacian determines ranking scores, a factor graph whose sparse Jacobian yields a sparse normal matrix, a communication graph over which nodes cooperatively solve an over-determined linear system, a survey network described by nodes and measurement edges, or a graph supporting time-varying signals and spectral filtering. Across these variants, the common operation is the minimization of squared residuals or weighted squared residuals, while the graph determines locality, sparsity, propagation, or gluing constraints (Hirani et al., 2010, Dong et al., 2019, Yang et al., 2018, Yan et al., 2024, Agbachi, 2018, Glass, 6 Mar 2026).

1. Graph structure as the organizing principle

A recurring feature of least-squares graph formulations is that the unknowns and the data inhabit different parts of a graph-derived object. In the ranking formulation on an oriented graph, the data are a 1-cochain $\omega \in C^1(G)$ on edges, while the unknown ranking is a 0-cochain $\alpha \in C^0(G)$ on vertices, and the model is

$\partial_1^T \alpha \simeq \omega,$

with least-squares fit

$\min_a \|\omega-\partial_1^T a\|_2.$

In factor-graph optimization, by contrast, the graph indexes residual factors over subsets of variables on a manifold, and the objective is

$x^{*}=\argmin_x \sum_i \|h_i(x)\|^2,$

which linearizes to

$$\Delta x^{*}=\argmin_{\Delta x}\|J\Delta x+b\|^2, \qquad J^T J\,\Delta x^{*}=J^T b.$$

In distributed solvers for network linear equations, each node holds one row $h_i^\top y=z_i$ of an over-determined system and updates a local state by combining consensus terms with local least-squares gradients. In online graph signal estimation, the graph supplies a Laplacian eigenspace and the least-squares signal update is driven by the residual observed on sampled nodes (Hirani et al., 2010, Dong et al., 2019, Liu et al., 2018, Yan et al., 2024).

Setting	Graph object	Canonical least-squares form
Ranking from pairwise comparisons	Weighted, oriented graph or comparison multigraph	$\min_a \|\omega-\partial_1^T a\|_2$; $Lq=s$
Robotics and vision	Factor graph	$\argmin_x \sum_i \|h_i(x)\|^2$
Distributed linear equations	Communication graph	$\alpha \in C^0(G)$0 with Laplacian-coupled dynamics
Online graph signals	Graph with Laplacian $\alpha \in C^0(G)$1	$\alpha \in C^0(G)$2

This range of meanings matters because the word graph does not play a single role. In some formulations it is the domain on which data live; in others it is the sparsity pattern of a nonlinear least-squares problem; in yet others it is the communication substrate of a distributed algorithm or the combinatorial scaffold for topological gluing. A plausible implication is that the phrase denotes a methodological family unified by residual minimization and graph-induced structure rather than by one canonical matrix equation.

2. Ranking, Laplacians, and Hodge decomposition

The most classical graph-theoretic least-squares formulation in the supplied literature concerns ranking from pairwise comparisons. On a weighted, oriented graph treated as an oriented abstract simplicial 1-complex, the unknown ranking is a vertex potential and the observed comparisons are edge values. Exact consistency would require the edge data to be a gradient field; when signed sums around cycles do not vanish, one instead computes the least-squares projection onto $\alpha \in C^0(G)$3. The normal equations are

$\alpha \in C^0(G)$4

where

$\alpha \in C^0(G)$5

is the graph Laplacian. This formulation reveals that ranking is a graph Laplacian problem and connects naturally to spectral graph theory, Hodge decomposition, KKT systems, and Betti numbers (Hirani et al., 2010).

The same idea appears in generalized tournaments with missing and multiple paired comparisons. There, one has objects $\alpha \in C^0(G)$6, an additive paired comparison matrix $\alpha \in C^0(G)$7, a symmetric matches matrix $\alpha \in C^0(G)$8, and a score vector

$\alpha \in C^0(G)$9

The least-squares rating is defined by

$\partial_1^T \alpha \simeq \omega,$0

with normalization

$\partial_1^T \alpha \simeq \omega,$1

The first-order conditions yield the Laplacian system

$\partial_1^T \alpha \simeq \omega,$2

so that

$\partial_1^T \alpha \simeq \omega,$3

For a connected comparison multigraph,

$\partial_1^T \alpha \simeq \omega,$4

The least-squares rating is unique if and only if the comparison multigraph is connected; equivalently, $\partial_1^T \alpha \simeq \omega,$5 is not block diagonal, the graph is connected, and the second-smallest Laplacian eigenvalue is positive (Csató, 2015).

A notable contribution of the generalized-tournament formulation is the iterative decomposition

$\partial_1^T \alpha \simeq \omega,$6

with

$\partial_1^T \alpha \simeq \omega,$7

where $\partial_1^T \alpha \simeq \omega,$8. This gives a graph interpretation in which direct scores are corrected by one-step, two-step, and higher-order indirect effects. To write the decomposition cleanly, the paper introduces a balanced comparison multigraph by attaching loops so that node $\partial_1^T \alpha \simeq \omega,$9 gets $\min_a \|\omega-\partial_1^T a\|_2.$0 loops. The iteration converges if the comparison graph is connected and not regular bipartite (Csató, 2015).

The Hodge-theoretic extension refines the residual further. When triangles are present, edge data admit the decomposition

$\min_a \|\omega-\partial_1^T a\|_2.$1

where $\min_a \|\omega-\partial_1^T a\|_2.$2 is the global ranking signal, $\min_a \|\omega-\partial_1^T a\|_2.$3 captures local triangular inconsistency, and $\min_a \|\omega-\partial_1^T a\|_2.$4 captures topological inconsistency. In this view, least-squares graph ranking is not only a ranking method but also an inconsistency analysis on a simplicial complex, with

$\min_a \|\omega-\partial_1^T a\|_2.$5

and $\min_a \|\omega-\partial_1^T a\|_2.$6 (Hirani et al., 2010).

3. Factor graphs, sparse nonlinear least squares, and robust graph-SLAM

In robotics and computer vision, least-squares graphs are most often factor graphs. miniSAM formulates the core problem on a manifold $\min_a \|\omega-\partial_1^T a\|_2.$7 as

$\min_a \|\omega-\partial_1^T a\|_2.$8

with Mahalanobis norm

$\min_a \|\omega-\partial_1^T a\|_2.$9

and after factorizing $x^{*}=\argmin_x \sum_i \|h_i(x)\|^2,$0 and dropping the robust loss for simplicity,

$x^{*}=\argmin_x \sum_i \|h_i(x)\|^2,$1

Linearization around $x^{*}=\argmin_x \sum_i \|h_i(x)\|^2,$2 gives

$x^{*}=\argmin_x \sum_i \|h_i(x)\|^2,$3

and, in stacked form,

$x^{*}=\argmin_x \sum_i \|h_i(x)\|^2,$4

The block sparsity of $x^{*}=\argmin_x \sum_i \|h_i(x)\|^2,$5 and $x^{*}=\argmin_x \sum_i \|h_i(x)\|^2,$6 follows directly from factor-graph locality: each factor touches only a small subset of variables (Dong et al., 2019).

The probabilistic interpretation is explicit. A factor graph is written as

$x^{*}=\argmin_x \sum_i \|h_i(x)\|^2,$7

and if each factor is Gaussian over its residual,

$x^{*}=\argmin_x \sum_i \|h_i(x)\|^2,$8

then MAP estimation is equivalent to minimizing a sum of squared residuals. miniSAM operationalizes this equivalence through a FactorGraph container, a Variables container, built-in PriorFactor and BetweenFactor, non-linear optimizers such as Levenberg–Marquardt, custom Factor subclasses with error() and optional jacobians(), and NumericalFactor for finite-difference derivatives. It includes built-in support for vector spaces and the Lie groups $x^{*}=\argmin_x \sum_i \|h_i(x)\|^2,$9, $$\Delta x^{*}=\argmin_{\Delta x}\|J\Delta x+b\|^2, \qquad J^T J\,\Delta x^{*}=J^T b.$$0, $$\Delta x^{*}=\argmin_{\Delta x}\|J\Delta x+b\|^2, \qquad J^T J\,\Delta x^{*}=J^T b.$$1, $$\Delta x^{*}=\argmin_{\Delta x}\|J\Delta x+b\|^2, \qquad J^T J\,\Delta x^{*}=J^T b.$$2, and $$\Delta x^{*}=\argmin_{\Delta x}\|J\Delta x+b\|^2, \qquad J^T J\,\Delta x^{*}=J^T b.$$3, and supports sparse backends including Eigen’s simplicial LDLT, CHOLMOD, and a CUDA-enabled cuSOLVER Cholesky solver. In the reported benchmarks, CHOLMOD gives the best CPU performance among the tested options, whereas CUDA cuSOLVER is not competitive on small problems and suffers from a one-time launch overhead of about 350 ms; GPU sparse solving is therefore only attractive for sufficiently large problems in the current implementation (Dong et al., 2019).

Graph-based SLAM introduces a second theme: robustification of least squares without abandoning factor-graph machinery. AEROS starts from the standard graph-SLAM objective and replaces the fixed robust kernel or per-loop switch with a single continuous latent parameter $$\Delta x^{*}=\argmin_{\Delta x}\|J\Delta x+b\|^2, \qquad J^T J\,\Delta x^{*}=J^T b.$$4 based on Barron’s general adaptive robust loss,

$$\Delta x^{*}=\argmin_{\Delta x}\|J\Delta x+b\|^2, \qquad J^T J\,\Delta x^{*}=J^T b.$$5

where $$\Delta x^{*}=\argmin_{\Delta x}\|J\Delta x+b\|^2, \qquad J^T J\,\Delta x^{*}=J^T b.$$6 gives ordinary quadratic/L2 loss, $$\Delta x^{*}=\argmin_{\Delta x}\|J\Delta x+b\|^2, \qquad J^T J\,\Delta x^{*}=J^T b.$$7 gives Cauchy-like loss, $$\Delta x^{*}=\argmin_{\Delta x}\|J\Delta x+b\|^2, \qquad J^T J\,\Delta x^{*}=J^T b.$$8 gives Welsh/Leclerc-type behavior, and $$\Delta x^{*}=\argmin_{\Delta x}\|J\Delta x+b\|^2, \qquad J^T J\,\Delta x^{*}=J^T b.$$9 gives pseudo-Huber / L1-L2 behavior. Using Black-Rangarajan duality, the robust penalty is written as

$h_i^\top y=z_i$0

and the joint optimization becomes

$h_i^\top y=z_i$1

Because the outlier-process term is nonnegative, it can be represented as a squared residual, making the whole construction compatible with Gaussian-factor solvers such as iSAM and iSAM2 (Ramezani et al., 2021).

The experimental protocol for AEROS uses synthetic 2D datasets Manhattan3500, CSAIL, INTEL, City10000, a synthetic 3D dataset Sphere2500, and the real 3D LiDAR dataset Newer College. The synthetic tests use 10 Monte Carlo trials per dataset with false loop closures injected at outlier ratios from 10% to 50%. The metric is Absolute Translation Error (ATE), and for Newer College the evaluation additionally uses Relative Translation Error (RTE) with Umeyama alignment. The reported pattern is that AEROS is stable across increasing outlier rates, is competitive with SC, DCS, GNC, and Geman-McClure (GM), and is especially strong on Newer College because it can softly absorb partially correct loop closures rather than forcing a binary inlier–outlier decision (Ramezani et al., 2021).

4. Distributed least squares over communication graphs

A different meaning of least-squares graph arises when the graph is a communication network. In the continuous-time formulation for network linear equations, each node knows only its local equation

$h_i^\top y=z_i$2

and the global objective is

$h_i^\top y=z_i$3

Each node maintains $h_i^\top y=z_i$4 and evolves according to

$h_i^\top y=z_i$5

where the first term is consensus and the second is local descent. Under the step-size conditions

$h_i^\top y=z_i$6

a fixed, connected graph and $h_i^\top y=z_i$7 yield convergence to the unique least-squares solution

$h_i^\top y=z_i$8

The paper further gives explicit rates: for $h_i^\top y=z_i$9 the rate depends on the curvature ratio $\min_a \|\omega-\partial_1^T a\|_2$0, and for $\min_a \|\omega-\partial_1^T a\|_2$1 with $\min_a \|\omega-\partial_1^T a\|_2$2 the average-state error decays as $\min_a \|\omega-\partial_1^T a\|_2$3 (Liu et al., 2018).

The Arrow-Hurwicz-Uzawa formulation sharpens the graph dependence. For the constrained problem

$\min_a \|\omega-\partial_1^T a\|_2$4

with

$\min_a \|\omega-\partial_1^T a\|_2$5

the continuous-time primal-dual flow is

$\min_a \|\omega-\partial_1^T a\|_2$6

Here connectivity alone is not sufficient. The paper proves a necessary-and-sufficient graph/data compatibility criterion: for every Laplacian eigenvector $\min_a \|\omega-\partial_1^T a\|_2$7, the set of row vectors $\min_a \|\omega-\partial_1^T a\|_2$8 on its support must span $\min_a \|\omega-\partial_1^T a\|_2$9; otherwise undamped oscillatory modes may remain. This condition explains why path and ring graphs can be favorable, while star graphs and complete graphs satisfy

$Lq=s$0

so for $Lq=s$1 the sufficient condition generally fails (Liu et al., 2017).

The discrete-time Euler discretization preserves the same structure:

$Lq=s$2

with stability threshold

$Lq=s$3

For $Lq=s$4, node states converge exponentially to $Lq=s$5; for $Lq=s$6, divergence can occur. The same paper also reports a distinctive switching-graph phenomenon: sufficiently fast switching can lead to approximate least square solvers even if all graphs in the switching signal fail to do so individually (Liu et al., 2017).

A discrete gradient-tracking alternative achieves exact exponential convergence on connected undirected graphs and strongly connected directed graphs. For the undirected case,

$Lq=s$7

$Lq=s$8

with $Lq=s$9. The paper proves the exact threshold

$\argmin_x \sum_i \|h_i(x)\|^2$0

which is both necessary and sufficient for exponential convergence. For strongly connected digraphs the corresponding push-pull method converges for sufficiently small $\argmin_x \sum_i \|h_i(x)\|^2$1. A separate finite-time mechanism reconstructs the exact least-squares solution from a single node’s local trajectory by detecting singularity in a Hankel matrix of state differences; in the reported 4-node examples, the finite-time solver recovers the exact solution in 16 time steps, whereas the plain iterative method needs around 300 steps to approach the solution numerically (Yang et al., 2018).

5. Matrix construction, large datasets, survey networks, and online graph signals

Least-squares graph methods also have a strongly computational interpretation. In polynomial curve fitting, the least-squares objective

$\argmin_x \sum_i \|h_i(x)\|^2$2

leads to the normal-equation system

$\argmin_x \sum_i \|h_i(x)\|^2$3

where $\argmin_x \sum_i \|h_i(x)\|^2$4 is built from sums of powers of the input $\argmin_x \sum_i \|h_i(x)\|^2$5-values and $\argmin_x \sum_i \|h_i(x)\|^2$6 from sums of $\argmin_x \sum_i \|h_i(x)\|^2$7. The paper’s contribution is a matricized formulation intended to expose parallelism: the costly part for large $\argmin_x \sum_i \|h_i(x)\|^2$8 is not the final $\argmin_x \sum_i \|h_i(x)\|^2$9 solve but the accumulation of the moment sums needed to build $\alpha \in C^0(G)$00 and $\alpha \in C^0(G)$01. The implementation uses CUDA on an NVIDIA Quadro 4000 GPU with 256 cores, reports speed-ups on the order of about $\alpha \in C^0(G)$02 for datasets with thousands of points when compared with sequential execution on a conventional multi-core processor, and solves the system by Gaussian elimination rather than explicit inversion (Dasgupta, 2015).

That paper also makes the numerical-robustness issue explicit by comparing with MATLAB’s polyfit(), which uses a Vandermonde matrix $\alpha \in C^0(G)$03 and QR factorization,

$\alpha \in C^0(G)$04

so that

$\alpha \in C^0(G)$05

The reported coefficients closely match polyfit() for linear, quadratic, and cubic fits, with correlation values around $\alpha \in C^0(G)$06 to $\alpha \in C^0(G)$07. For the sample cubic case, the summed squared errors are

$\alpha \in C^0(G)$08

and the paper therefore claims that the generated coefficients produce a best-fit curve according to the least-squares criterion (Dasgupta, 2015).

In Geomatics Engineering, least-squares graphs appear as survey networks. The adjustment model is

$\alpha \in C^0(G)$09

and after introducing approximate values $\alpha \in C^0(G)$10,

$\alpha \in C^0(G)$11

Weighted least squares then minimizes

$\alpha \in C^0(G)$12

yielding the normal equations

$\alpha \in C^0(G)$13

and solution

$\alpha \in C^0(G)$14

with covariance

$\alpha \in C^0(G)$15

The network is viewed as a directed graph in which nodes are stations and edges are observed directions or measurements. The paper emphasizes cycle processing through Breadth First Search (BFS) and Depth First Search (DFS), preferring DFS because survey routes naturally follow depth-like progression and DFS has lower memory requirements. Its computational claim is that the manual stages of Collation, Classification, Data input, and Computation can be replaced by a frame-based and object-oriented workflow containing operations such as FormEquations and SolveForX (Agbachi, 2018).

Online graph signal estimation supplies a further computational extension. LMS-GNN considers noisy, partially observed graph signals

$\alpha \in C^0(G)$16

and uses the mean-squared objective

$\alpha \in C^0(G)$17

With residual

$\alpha \in C^0(G)$18

the LMS-style update is

$\alpha \in C^0(G)$19

Its deep version uses

$\alpha \in C^0(G)$20

On hourly temperatures from $\alpha \in C^0(G)$21 U.S. weather stations over $\alpha \in C^0(G)$22 time steps, with the first 24 used for training, an 8-nearest neighbors graph, and Gaussian noise variance $\alpha \in C^0(G)$23, LMS-GNN achieves the best reported Spatial MSE and Spectral MAE among GLMS, GNLMS, GCN, and STGCN. For example, the Spatial MSE at $\alpha \in C^0(G)$24 is $\alpha \in C^0(G)$25 for LMS-GNN versus $\alpha \in C^0(G)$26, $\alpha \in C^0(G)$27, $\alpha \in C^0(G)$28, and $\alpha \in C^0(G)$29 for GLMS, GNLMS, GCN, and STGCN, respectively (Yan et al., 2024).

6. Random-walk, geometric, and homotopy-theoretic extensions

A more recent reinterpretation recasts least squares as a cumulative graphical process. For equally spaced samples $\alpha \in C^0(G)$30, define the mean-centered cumulative sum

$\alpha \in C^0(G)$31

which produces a pinned data walk satisfying $\alpha \in C^0(G)$32. The signed area under the walk is

$\alpha \in C^0(G)$33

For uniformly spaced design points

$\alpha \in C^0(G)$34

the paper proves that the slope which annuls the net signed area under the residual data walk is exactly the ordinary least-squares slope. In its notation,

$\alpha \in C^0(G)$35

and this coincides with the conventional LLS formula for equally spaced samples. The result is exact for arbitrary distributions of steps; for non-uniform sampling it is approximate rather than exact (Kostinski et al., 26 Mar 2025).

A separate geometric line of work, while not graph-theoretic in the combinatorial sense, gives an adjacent interpretation of least-squares line fitting through covariance geometry. When both variables have all equal uncorrelated errors, the least-squares fit minimizes the sum of squared perpendicular distances to a line, and the best-fit direction is the principal component direction of maximal variance of the data cloud. In angle form, the fit satisfies

$\alpha \in C^0(G)$36

with closed-form uncertainty formulas for $\alpha \in C^0(G)$37, $\alpha \in C^0(G)$38, slope, and intercept. This perspective is important because it isolates a regime in which least-squares line fitting reduces to a covariance-matrix eigenproblem (Petrolini, 2011).

The most abstract extension in the supplied literature is homotopy-theoretic least squares regression. There, weighted finite datasets form a category $\alpha \in C^0(G)$39, each object carries a Koszul complex built from the least-squares normal equations, and the zeroth homology

$\alpha \in C^0(G)$40

recovers the coordinate ring of the least-squares solution locus. To compare local least-squares solutions on overlapping subsets, the paper linearizes near chosen solutions by passing to

$\alpha \in C^0(G)$41

and restores functoriality by translation maps

$\alpha \in C^0(G)$42

Evaluating the resulting presheaf on a cover yields a Čech-Koszul bicomplex in which degree-0 cocycles are local least-squares solutions, degree-1 terms are homotopies on overlaps, and higher-degree terms are higher coherences. The paper explicitly states that the resulting “least-squares graph” has local LS solutions as vertices and edges or higher simplices encoding discrepancies between them via chain homotopies (Glass, 6 Mar 2026).

This broad set of formulations shows that least-squares graph methods range from classical Laplacian ranking and sparse factor-graph optimization to distributed consensus solvers, adaptive graph filters, random-walk reformulations, and homotopical gluing constructions. A plausible implication is that the enduring role of graph structure is not merely representational: it determines what counts as locality, which residuals may be coupled, how inconsistency propagates, and which computational methods—QR, Gaussian elimination, sparse Cholesky, CUDA parallelism, DFS/BFS cycle processing, or distributed gradient tracking—are natural for the least-squares problem at hand.