Case Study 2: Building a Least-Squares Curve Fitter

"Fit the model you can defend, with the numerics you can trust."

Executive Summary

Fitting a polynomial to data is one of the most common jobs in a scientist's day: you have measurements and you want the smooth curve that best explains them. This study builds a small, reusable polynomial-fitting module from the pieces this chapter gave you — the transpose and matmul intrinsics to assemble the normal equations, and dgesv to solve them — wraps it behind a clean polyfit(x, y, deg, coef) interface, and then verifies it the only honest way: by handing it data sampled from a polynomial we already know and checking that it recovers the coefficients exactly. Finally it delivers the warning that separates a toolbox routine from a trap: the normal-equations approach we build is the readable one and the numerically weakest one, because it squares the condition number of the problem — so we end by showing when to reach past dgesv for the least-squares driver dgels or the SVD instead. You will finish with a working fitter and, more importantly, with the judgment to know its limits.

Skills applied: designing a reusable numerical routine with a clean interface; assembling the Vandermonde matrix and normal equations with transpose/matmul (§21.1); solving with dgesv and checking info (§21.3); the leading dimension for allocatable matrices (§21.3); the conditioning cost of the normal equations (§21.4 and Chapter 20); choosing dgels/SVD for robustness (§21.4).

Background

You are given $m$ data points $(x_i, y_i)$ and want the degree-$d$ polynomial $p(x) = c_1 + c_2 x + \dots + c_{d+1} x^{d}$ whose values come closest, in the least-squares sense, to the $y_i$. Stack the "design" (Vandermonde) matrix $A$, whose $i$-th row is $[1,\ x_i,\ x_i^2,\ \dots,\ x_i^d]$, so that $A\mathbf{c} \approx \mathbf{y}$ is an overdetermined system — more equations ($m$) than unknowns ($d+1$). The least-squares solution is the $\mathbf{c}$ minimizing $\lVert A\mathbf{c} - \mathbf{y}\rVert^2$, and the classic route to it is the normal equations

$$ A^{\mathsf{T}} A\,\mathbf{c} = A^{\mathsf{T}}\mathbf{y}, $$

a square $(d{+}1)\times(d{+}1)$ system you can hand straight to dgesv. Everything you need is in this chapter: transpose and matmul build $A^{\mathsf{T}}A$ and $A^{\mathsf{T}}\mathbf{y}$, and dgesv solves.

Phase 1 — Design the Interface

Before writing a line of arithmetic, decide what the routine promises. A good fitter hides all the linear algebra and exposes only the data, the degree, and the answer:

!  polyfit(x, y, deg, coef, info)
!    x, y  : the data points (length m)
!    deg   : polynomial degree; coef has deg+1 entries
!    coef  : OUTPUT, the coefficients [c0, c1, ..., c_deg]
!    info  : OUTPUT, 0 on success (passed straight through from dgesv)

Passing info out rather than stopping inside is a deliberate design choice: a library routine reports failure and lets the caller decide what to do, exactly as dgesv reports to us. The routine allocates its own workspace from deg, so the caller never sees a Vandermonde matrix or a pivot array.

Phase 2 — Assemble the Normal Equations

The Vandermonde matrix has one column per power of $x$: column $j$ is $x^{\,j-1}$, which a whole-array exponentiation writes in one line. Then $A^{\mathsf{T}}A$ and $A^{\mathsf{T}}\mathbf{y}$ are two matmuls with a transpose — the §21.1 intrinsics doing exactly the linear algebra they are for:

do j = 1, ncoef
  vand(:, j) = x**(j - 1)                 ! columns: 1, x, x^2, ...  (x**0 = 1)
end do
ata = matmul(transpose(vand), vand)       ! (d+1)x(d+1) normal-equations matrix
aty = matmul(transpose(vand), y)          ! (d+1) right-hand side

Note the shapes line up because transpose(vand) is $(d{+}1)\times m$ and vand is $m\times(d{+}1)$, so their product is the square $(d{+}1)\times(d{+}1)$ that dgesv wants. This is why matrix–matrix intrinsics earn their place: the alternative is a double loop of accumulations you would have to index by hand.

Phase 3 — Build the Module and Solve

Here is the complete, self-contained fitter. The leading dimension passed to dgesv is ncoef, the declared first dimension of the square ata — the point from §21.3 in a routine with allocatable arrays:

module polyfit_mod
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
contains
  subroutine polyfit(x, y, deg, coef, info)
    real(dp), intent(in)  :: x(:), y(:)
    integer,  intent(in)  :: deg
    real(dp), intent(out) :: coef(deg+1)
    integer,  intent(out) :: info
    integer               :: ncoef, j
    real(dp), allocatable :: vand(:,:), ata(:,:), aty(:)
    integer,  allocatable :: ipiv(:)

    ncoef = deg + 1
    allocate(vand(size(x), ncoef), ata(ncoef, ncoef), aty(ncoef), ipiv(ncoef))
    do j = 1, ncoef
      vand(:, j) = x**(j - 1)
    end do
    ata = matmul(transpose(vand), vand)
    aty = matmul(transpose(vand), y)
    call dgesv(ncoef, 1, ata, ncoef, ipiv, aty, ncoef, info)   ! aty <- coefficients
    coef = aty
  end subroutine polyfit
end module polyfit_mod

program fit_demo
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use polyfit_mod
  implicit none
  real(dp) :: x(4) = [ -1.0_dp, 0.0_dp, 1.0_dp, 2.0_dp ]
  real(dp) :: y(4) = [  0.0_dp, 2.0_dp, 6.0_dp, 12.0_dp ]    ! sampled from 2 + 3x + x^2
  real(dp) :: coef(3)
  integer  :: info

  call polyfit(x, y, 2, coef, info)
  if (info /= 0) then
    print '(a, i0)', 'fit failed, info = ', info
    error stop 1
  end if
  print '(a)', 'fitted coefficients [c0, c1, c2]:'
  print '(f8.4)', coef
end program fit_demo

$ gfortran -std=f2018 -Wall polyfit.f90 -o fit -llapack -lblas && ./fit
fitted coefficients [c0, c1, c2]:
  2.0000
  3.0000
  1.0000

(The module precedes the program in the single source file, so gfortran compiles them in the right order; the driver uses polyfit_mod.)

Phase 4 — Verify Against a Known Answer

A fitter that returns numbers is not a fitter you trust; a fitter that recovers a polynomial you chose is. The data was sampled exactly from $p(x) = 2 + 3x + x^2$: at $x = -1, 0, 1, 2$ that gives $y = 0, 2, 6, 12$ (check $x=2$: $2 + 6 + 4 = 12$). Because the points lie exactly on a quadratic, the least-squares fit is exact, and the coefficients must come back as $(2, 3, 1)$ — which they do. Work the normal equations by hand to see the machinery is right. The Vandermonde columns give the sums

$$ A^{\mathsf{T}}A = \begin{bmatrix} 4 & 2 & 6 \\ 2 & 6 & 8 \\ 6 & 8 & 18 \end{bmatrix}, \qquad A^{\mathsf{T}}\mathbf{y} = \begin{bmatrix} 20 \\ 30 \\ 54 \end{bmatrix}, $$

where, for instance, the top-left $4$ is $\sum 1 = 4$, the $6$ in the corner is $\sum x^2 = 1+0+1+4 = 6$, and the last right-hand entry is $\sum x^2 y = 0 + 0 + 6 + 48 = 54$. Substituting $\mathbf{c} = (2, 3, 1)$ into the first row: $4\cdot2 + 2\cdot3 + 6\cdot1 = 8 + 6 + 6 = 20$ — matches. The routine reproduces the polynomial it was fed, so the assembly and the solve are both correct. This "recover a known model" check is the unit test every fitting routine should ship with.

Phase 5 — The Catch: Normal Equations Square the Conditioning

The routine is correct, readable, and — for high degrees or clustered data — the least accurate way to fit. The reason is exactly the conditioning idea from Chapter 20, sharpened by a fact about the normal equations: forming $A^{\mathsf{T}}A$ squares the condition number. If the design matrix $A$ has condition number $\kappa(A)$, then $\kappa(A^{\mathsf{T}}A) \approx \kappa(A)^2$. A Vandermonde matrix on a wide range of $x$ values is already ill-conditioned — its columns $1, x, x^2, \dots$ grow nearly parallel — so a degree-8 fit whose $A$ has $\kappa(A) \approx 10^{6}$ hands dgesv a normal-equations matrix with $\kappa \approx 10^{12}$, and in double precision (about 16 digits) you have thrown away twelve of them before the solve begins. The fit still returns numbers; they are just increasingly meaningless as the degree climbs.

The professional fix is to not form the normal equations at all. LAPACK's dgels solves the least-squares problem directly from $A$ by a QR factorization, and the SVD routine dgesvd (§21.4) does it even more robustly and tells you, through the singular values, exactly how rank-deficient the fit is. Both work on the design matrix $A$ itself, so they see $\kappa(A)$, not $\kappa(A)^2$ — the difference between six trustworthy digits and none. Our dgesv-on-normal-equations fitter is the right tool for low degrees and well-spread data, where it is simplest and perfectly accurate; the moment the degree grows or the fit looks unstable, you swap the solver — the same "call the right specialist" judgment §21.6 made for sparse matrices, now for least squares.

Approach What it factors Conditioning it sees Use when
Normal equations + dgesv (this build) $A^{\mathsf{T}}A$ $\kappa(A)^2$ low degree, well-spread data; simplicity matters
QR + dgels $A$ directly $\kappa(A)$ the default robust choice for least squares
SVD + dgesvd $A$ directly $\kappa(A)$, and reports rank ill-conditioned or rank-deficient fits; you need diagnostics

Discussion Questions

  1. The build allocates vand, ata, aty, and ipiv from deg on every call. For a routine called thousands of times in a loop, what would you change, and what does intent(out) on coef already guarantee about aliasing that makes the change safe?
  2. We verified the fitter by recovering a polynomial from exact samples. Why is that a stronger test than fitting real noisy data and eyeballing the curve — and what property of the least-squares solution makes the exact-recovery test pass to the last digit?
  3. dgels avoids squaring the condition number but needs a workspace query (like dsyev) and returns the solution in the first n entries of the right-hand-side array. Sketch how the polyfit interface would change if you swapped dgesv for dgels, and what would not change for the caller.

Your Turn: Extensions

  • Option A. Add a residual and goodness-of-fit report: after solving, evaluate the fitted polynomial at each x (Horner's rule, or a second Vandermonde matmul), and print the maximum residual $\max_i |p(x_i) - y_i|$. On the exact-quadratic data it should be at the rounding floor; verify that.
  • Option B. Push the degree until it breaks. Fit degree 2, 4, 6, 8 to points sampled from a known high-degree polynomial on $x \in [0, 10]$, and watch the recovered coefficients drift from the true ones as the normal-equations conditioning degrades. Where does double precision give out?
  • Option C. Rebuild polyfit on dgels instead of the normal equations (workspace query and all), and repeat Option B. Show that the QR-based fit holds accuracy several degrees further — the payoff of seeing $\kappa(A)$ instead of $\kappa(A)^2$.

Key Takeaways

  • A least-squares polynomial fit is transpose + matmul to build the normal equations and one dgesv to solve them — the chapter's intrinsics and driver composed into a genuine, reusable tool.
  • Wrap the linear algebra behind a clean interface (polyfit(x, y, deg, coef, info)) and pass info out: a library routine reports failure and lets the caller decide, the way dgesv reports to you.
  • Verify a numerical routine by recovering a known answer — sample a polynomial you chose and confirm the fitter returns its coefficients. That is the unit test, not a plot you squint at.
  • The normal equations are the readable route and the numerically weakest one: they square the condition number (Chapter 20). For high degree or ill-conditioned data, solve the least-squares problem directly with dgels (QR) or dgesvd (SVD) — see $\kappa(A)$, not $\kappa(A)^2$.