---
title: Least-Squares Graph Methods
url: https://www.emergentmind.com/topics/least-squares-graph
type: topic
---

# Least-Squares Graph Methods

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 [1011.1716][1909.00903][1810.00156][2401.15304][1804.05665][2603.05870].

## 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 [1011.1716][1909.00903][1808.04140][2401.15304].

| 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 | $\min_y \|z-Hy\|^2$ with Laplacian-coupled dynamics |
| Online graph signals | Graph with Laplacian $\mathbf L$ | $\mathbb E\|\boldsymbol y[t]-\mathbf D_{\mathcal S}\mathbf B_N\hat{\boldsymbol x}[t]\|_2^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 $\operatorname{im}(\partial_1^T)$. The normal equations are
$$
\partial_1\partial_1^T a=\partial_1\omega,
$$
where
$$
\Delta_0=\partial_1\partial_1^T
$$
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 [1011.1716].

The same idea appears in generalized tournaments with missing and multiple paired comparisons. There, one has objects $N=\{X_1,\dots,X_n\}$, an additive paired comparison matrix $A=(a_{ij})$, a symmetric matches matrix $M=(m_{ij})$, and a score vector
$$
s=Ae.
$$
The least-squares rating is defined by
$$
\min_{q\in\mathbb R^n}\sum_{X_i,X_j\in N} m_{ij}(h_{ij}-q_i+q_j)^2,
$$
with normalization
$$
e^Tq=0.
$$
The first-order conditions yield the Laplacian system
$$
Lq=s,
$$
so that
$$
q=L^+s, \qquad e^Tq=0.
$$
For a connected comparison multigraph,
$$
q=L^+s=\left[L+\frac{1}{n}J\right]^{-1}s.
$$
The least-squares rating is unique if and only if the comparison multigraph is connected; equivalently, $M$ is not block diagonal, the graph is connected, and the second-smallest Laplacian eigenvalue is positive [1508.06778].

A notable contribution of the generalized-tournament formulation is the iterative decomposition
$$
q=\frac{1}{d}\sum_{k=0}^{\infty} P^k s,
$$
with
$$
q^{(0)}=\frac{1}{d}s,\qquad
q^{(k)}=q^{(k-1)}+\frac{1}{d}P^k s,
\qquad
q=\lim_{k\to\infty}q^{(k)},
$$
where $d=\max_i d_i$. 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 $X_i$ gets $d-d_i$ loops. The iteration converges if the comparison graph is connected and not regular bipartite [1508.06778].

The Hodge-theoretic extension refines the residual further. When triangles are present, edge data admit the decomposition
$$
\omega=\partial_1^T\alpha+\partial_2\beta+h,
$$
where $\partial_1^T\alpha$ is the global ranking signal, $\partial_2\beta$ captures local triangular inconsistency, and $h\in\ker\Delta_1$ 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
$$
C^1(G)=\operatorname{im}(\partial_1^T)\oplus \operatorname{im}(\partial_2)\oplus \ker(\Delta_1),
$$
and $\dim\ker\Delta_p=\beta_p$ [1011.1716].

## 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 $\mathcal M$ as
$$
x^{*}=\argmin_x \sum_i \rho_i\big(\|f_i(x)\|_{\Sigma_i}^2\big),
$$
with Mahalanobis norm
$$
\|v\|^2_{\Sigma}=v^T\Sigma^{-1}v,
$$
and after factorizing $\Sigma^{-1}=R^TR$ and dropping the robust loss for simplicity,
$$
x^{*}=\argmin_x \sum_i \|h_i(x)\|^2.
$$
Linearization around $x_0$ gives
$$
\Delta x^{*}=\argmin_{\Delta x}\sum_i \|J_i\Delta x+h_i(x_0)\|^2,
\qquad
x^{*}=x_0\oplus \Delta x^{*},
$$
and, in stacked form,
$$
J^TJ\,\Delta x^{*}=J^Tb.
$$
The block sparsity of $J$ and $J^TJ$ follows directly from factor-graph locality: each factor touches only a small subset of variables [1909.00903].

The probabilistic interpretation is explicit. A factor graph is written as
$$
p(x)\propto \prod_i p_i(x_i),
$$
and if each factor is Gaussian over its residual,
$$
p_i(x_i)\propto \exp\!\left(-\frac12\|f_i(x_i)\|_{\Sigma_i}^2\right),
$$
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 $SO(2)$, $SE(2)$, $SO(3)$, $SE(3)$, and $Sim(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 [1909.00903].

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** $\alpha$ based on Barron’s general adaptive robust loss,
$$
\rho(\nu,\alpha,c)=\frac{\alpha-2}{\alpha}
\left(
\left(\frac{(\nu/c)^2}{\alpha-2}+1\right)^{\alpha/2}-1
\right),
$$
where $\alpha\to 2$ gives ordinary quadratic/L2 loss, $\alpha\to 0$ gives Cauchy-like loss, $\alpha\to -\infty$ gives Welsh/Leclerc-type behavior, and $\alpha=1$ gives pseudo-Huber / L1-L2 behavior. Using Black-Rangarajan duality, the robust penalty is written as
$$
\rho(\nu_j)=\omega_j\nu_j^2+\Psi(\omega_j),
$$
and the joint optimization becomes
$$
\argmin_{\mathcal X,\alpha}
\sum_i \nu_i^2+\sum_j\big(\omega_j\nu_j^2+\Psi(\omega_j,\alpha)\big).
$$
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 [2110.02018].

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 [2110.02018].

## 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,
$$
and the global objective is
$$
\min_{y\in\mathbb R^m}\|z-Hy\|^2
=
\min_{y\in\mathbb R^m}\sum_{i=1}^N |h_i^\top y-z_i|^2.
$$
Each node maintains $x_i(t)\in\mathbb R^m$ and evolves according to
$$
\dot{x}_i(t)
=
K\sum_{j\in\mathcal N_i(t)}A_{ij}(t)\big(x_j(t)-x_i(t)\big)
-\frac{\alpha(t)}{2}\nabla f_i(x_i(t)),
$$
where the first term is consensus and the second is local descent. Under the step-size conditions
$$
\int_0^\infty \alpha(t)\,dt=\infty,\qquad
\lim_{t\to\infty}\alpha(t)=0,\qquad
\int_0^\infty \alpha^2(t)\,dt<\infty,
$$
a fixed, connected graph and $\operatorname{rank}(H)=m$ yield convergence to the unique least-squares solution
$$
y^\ast=(H^\top H)^{-1}H^\top z.
$$
The paper further gives explicit rates: for $\alpha(t)=O(1/t)$ the rate depends on the curvature ratio $\frac{H^\top H}{N}$, and for $\alpha(t)=O(1/t^\lambda)$ with $\lambda\in(0,1)$ the average-state error decays as $O(t^{-\lambda})$ [1808.04140].

The Arrow-Hurwicz-Uzawa formulation sharpens the graph dependence. For the constrained problem
$$
\min_x U(x)\quad \text{s.t.}\quad (L\otimes I_m)x=0,
$$
with
$$
U(x_1,\dots,x_N)=\sum_{i=1}^N |h_i^\top x_i-z_i|^2,
$$
the continuous-time primal-dual flow is
$$
\dot x=-(L\otimes I_m)v-\nabla U(x),\qquad
\dot v=(L\otimes I_m)x.
$$
Here connectivity alone is not sufficient. The paper proves a necessary-and-sufficient graph/data compatibility criterion: for every Laplacian eigenvector $\alpha$, the set of row vectors $\{h_i:\ i\in\mathcal I_\alpha\}$ on its support must span $\mathbb R^m$; otherwise undamped oscillatory modes may remain. This condition explains why path and ring graphs can be favorable, while star graphs and complete graphs satisfy
$$
\min_{\alpha\in\mathcal S_L}|\mathcal I_\alpha|=2,
$$
so for $m>2$ the sufficient condition generally fails [1701.03908].

The discrete-time Euler discretization preserves the same structure:
$$
x(k+1)=x(k)-\epsilon(L\otimes I_m)v(k)-\epsilon\nabla U(x(k)),
\qquad
v(k+1)=v(k)+\epsilon(L\otimes I_m)x(k),
$$
with stability threshold
$$
\epsilon^\ast=
\min_{\lambda\in\sigma^\ast(M)}
\left[-\frac{2\Re(\lambda)}{|\lambda|^2}\right].
$$
For $0<\epsilon<\epsilon^\ast$, node states converge exponentially to $y^\ast$; for $\epsilon>\epsilon^\ast$, 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 [1701.03908].

A discrete gradient-tracking alternative achieves exact exponential convergence on connected undirected graphs and strongly connected directed graphs. For the undirected case,
$$
x_i(t+1)=\sum_{j\in\mathcal N_i}W_{ij}x_j(t)-\alpha v_i(t),
$$
$$
v_i(t+1)=\sum_{j\in\mathcal N_i}W_{ij}v_j(t)
+\nabla f_i(x_i(t+1))-\nabla f_i(x_i(t)),
$$
with $v_i(0)=\nabla f_i(x_i(0))$. The paper proves the exact threshold
$$
\alpha<\bar\alpha,\qquad
\bar\alpha=
\frac{1}{2\lambda_{\max}\!\left(\big((I_N+W)^{-2}\otimes I_m\big)\tilde H\right)},
$$
which is both necessary and sufficient for exponential convergence. For strongly connected digraphs the corresponding push-pull method converges for sufficiently small $\alpha>0$. 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 [1810.00156].

## 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
$$
I=\sum_{i=1}^{n}(y_i-f(x_i))^2,
\qquad
f(x)=a_0+a_1x+\cdots+a_mx^m,
$$
leads to the normal-equation system
$$
AX=B,
\qquad
X=
\begin{bmatrix}
a_0\\ a_1\\ a_2\\ \vdots\\ a_m
\end{bmatrix},
$$
where $A$ is built from sums of powers of the input $x$-values and $B$ from sums of $x_i^k y_i$. The paper’s contribution is a **matricized** formulation intended to expose parallelism: the costly part for large $n$ is not the final $(m+1)\times(m+1)$ solve but the accumulation of the moment sums needed to build $A$ and $B$. The implementation uses CUDA on an **NVIDIA Quadro 4000 GPU with 256 cores**, reports speed-ups on the order of about **$100\times$** 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 [1512.08017].

That paper also makes the numerical-robustness issue explicit by comparing with MATLAB’s `polyfit()`, which uses a Vandermonde matrix $V$ and QR factorization,
$$
Vp=Y,\qquad V=QR,\qquad Rp=Q^TY,
$$
so that
$$
p=R^{-1}Q^TY.
$$
The reported coefficients closely match `polyfit()` for linear, quadratic, and cubic fits, with correlation values around **$0.9996$ to $0.9998$**. For the sample cubic case, the summed squared errors are
$$
\sum e_f^2=128.199937,\qquad
\sum e_p^2=129.651164,
$$
and the paper therefore claims that the generated coefficients produce a best-fit curve according to the least-squares criterion [1512.08017].

In Geomatics Engineering, least-squares graphs appear as survey networks. The adjustment model is
$$
Y+v=AX+C,
$$
and after introducing approximate values $X=X_0+x$,
$$
v=Ax-b.
$$
Weighted least squares then minimizes
$$
v^TWv=\min,
$$
yielding the normal equations
$$
A^TWAx=A^TWb
$$
and solution
$$
x=(A^TWA)^{-1}A^TWb,
$$
with covariance
$$
\sigma_x^2=\sigma_0^2(A^TWA)^{-1},
\qquad
\sigma_0^2=\frac{v^TWv}{n-m}.
$$
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` [1804.05665].

Online graph signal estimation supplies a further computational extension. LMS-GNN considers noisy, partially observed graph signals
$$
\boldsymbol y[t]=\mathbf D_{\mathcal S}\big(\boldsymbol x_g[t]+\boldsymbol w[t]\big),
$$
and uses the mean-squared objective
$$
J(\hat{\boldsymbol x}[t])=
\mathbb E\left\|
\boldsymbol y[t]-\mathbf D_{\mathcal S}\mathbf B_N\hat{\boldsymbol x}[t]
\right\|_2^2.
$$
With residual
$$
\boldsymbol e[t]=\mathbf D_{\mathcal S}\big(\boldsymbol y[t]-\hat{\boldsymbol x}[t]\big),
$$
the LMS-style update is
$$
\hat{\boldsymbol x}[t+1]
=
\hat{\boldsymbol x}[t]+\mu_{\text{lms}}\mathbf B_N\boldsymbol e[t].
$$
Its deep version uses
$$
\hat{\boldsymbol x}_{l+1}[t]
=
\sigma\!\left(
\hat{\boldsymbol x}_l[t]
+
\mathbf U\sum_{i=1}^{f}\mathbf\Theta_{i,l}\mathbf U^\mathsf T\boldsymbol e_l[t]
+\boldsymbol b_l[t]
\right).
$$
On hourly temperatures from **$N=197$** U.S. weather stations over **$T=95$** time steps, with the first **24** used for training, an **8-nearest neighbors** graph, and Gaussian noise variance $\text{VAR}\in\{0.1,0.5,1\}$, LMS-GNN achieves the best reported **Spatial MSE** and **Spectral MAE** among **GLMS, GNLMS, GCN, and STGCN**. For example, the Spatial MSE at $\text{VAR}=0.1$ is **$0.555$** for LMS-GNN versus **$2.112$**, **$1.470$**, **$7.559$**, and **$4.990$** for GLMS, GNLMS, GCN, and STGCN, respectively [2401.15304].

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

A more recent reinterpretation recasts least squares as a cumulative graphical process. For equally spaced samples $y_k$, define the mean-centered cumulative sum
$$
z_j=\sum_{k=1}^j (y_k-\bar y),
\qquad
z_0=0,
$$
which produces a pinned **data walk** satisfying $z_N=0$. The signed area under the walk is
$$
A(y)=-\sum_{j=1}^N z_j.
$$
For uniformly spaced design points
$$
x_k=\frac{k-1}{N-1},
$$
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=
-\frac{12}{N(N+1)}\sum_{k=1}^N z_k,
$$
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 [2503.20965].

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
$$
\tan(2\theta)=\frac{2\,\mathrm{CovXY}}{\mathrm{VarX}-\mathrm{VarY}},
$$
with closed-form uncertainty formulas for $\hat\theta$, $\hat c$, slope, and intercept. This perspective is important because it isolates a regime in which least-squares line fitting reduces to a covariance-matrix eigenproblem [1104.3132].

The most abstract extension in the supplied literature is homotopy-theoretic least squares regression. There, weighted finite datasets form a category \(Fin\), each object carries a Koszul complex built from the least-squares normal equations, and the zeroth homology
$$
H_0(K_\bullet)\cong \mathcal R^{\omega\mathscr D}/(\eta^1,\dots,\eta^n)
$$
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
$$
\mathcal R^{\omega\mathscr D}_{\mathfrak a}
=
\mathcal R^{\omega\mathscr D}/\mathfrak I_{\mathfrak a}^2,
$$
and restores functoriality by translation maps
$$
\tau_{\mathfrak a,\mathfrak b}(a,\omega)=\bigl(a-(\mathfrak a-\mathfrak b),\omega\bigr).
$$
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 [2603.05870].

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.

Source: https://www.emergentmind.com/topics/least-squares-graph