The Geometry of Least Squares

Math 123: Mathematical Aspects of Data Analysis — Fall 2026

Published

September 29, 2026

Least squares is orthogonal projection. That is, the fitted values $\widehat{\bfy}$ are the orthogonal projection of $\bfy$ onto the column space $\Ccal(\bfX)$, and the residual $\widehat{\bfr}=\bfy-\widehat{\bfy}$ is perpendicular to $\Ccal(\bfX)$.

A tiny regression

Consider the setting of fitting a line $y=\theta_1+\theta_2x$ to the three points $(1,1)$, $(2,0)$, $(3,3)$. The design matrix is $\bfX=[\mathbf{1}\ \bfx]\in\Rbb^{3\times2}$, so its columns are the all-ones vector $\mathbf{1}$ and the feature vector $\bfx=(1,2,3)$.

import numpy as np
import plotly.graph_objects as go

x = np.array([1.0, 2.0, 3.0])
X = np.column_stack([np.ones(3), x])
y = np.array([1.0, 0.0, 3.0])

theta_hat, *_ = np.linalg.lstsq(X, y, rcond=None)
y_hat = X @ theta_hat
r_hat = y - y_hat

print("theta_hat =", np.round(theta_hat, 4))
print("y_hat     =", np.round(y_hat, 4))
print("r_hat     =", np.round(r_hat, 4))
print("X^T r_hat =", X.T @ r_hat)

theta_hat = [-0.6667  1.    ]
y_hat     = [0.3333 1.3333 2.3333]
r_hat     = [ 0.6667 -1.3333  0.6667]
X^T r_hat = [-6.66133815e-16 -1.77635684e-15]

The last line is the normal equations $\bfX^\top\widehat{\bfr}=\mathbf{0}$, up to rounding: the residual is orthogonal to both columns of $\bfX$, and hence to the whole plane $\Ccal(\bfX)$.

The picture

In Figure 1, each axis is one observation, not one feature. We are looking at $\Rbb^n$, the space where $\bfy$ lives. Drag to rotate, scroll to zoom, and click legend entries to hide pieces.

Plot the projection geometry

# The same muted palette as the other examples, which reads correctly against
# both the light and the dark version of this site.
INK = "#7c7c78"
GRID = "rgba(128, 128, 128, 0.28)"
BLUE, ORANGE, GREEN = "#2a78d6", "#eb6834", "#1baf7a"
ORIGIN = np.zeros(3)


def unit(v):
    return v / np.linalg.norm(v)


def vec(v):
    return "(" + ", ".join(f"{c:.3f}" for c in v) + ")"


def arrow(start, end, color, name, hover, width=7, head=0.25, showlegend=True, group=None):
    """A vector from start to end: a line capped by a cone, toggled together in the legend."""
    group = group or name
    d = unit(end - start)
    shaft = end - 0.5 * head * d  # tuck the end of the line inside the cone
    return [
        go.Scatter3d(x=[start[0], shaft[0]], y=[start[1], shaft[1]], z=[start[2], shaft[2]],
                     mode="lines", name=name, legendgroup=group, showlegend=showlegend,
                     line=dict(color=color, width=width),
                     hovertemplate=hover + "<extra></extra>"),
        go.Cone(x=[end[0]], y=[end[1]], z=[end[2]], u=[head * d[0]], v=[head * d[1]],
                w=[head * d[2]], anchor="tip", sizemode="raw", sizeref=1,
                colorscale=[[0, color], [1, color]], showscale=False,
                lighting=dict(ambient=0.9, diffuse=0.1, specular=0), hoverinfo="skip",
                legendgroup=group, showlegend=False),
    ]


def label(point, text, color, ax, ay):
    """A text label pinned to a 3D point, with a thin leader line and a border in the piece's color."""
    return dict(x=point[0], y=point[1], z=point[2], text=text, showarrow=True, arrowhead=0,
                arrowwidth=1, arrowcolor=INK, ax=ax, ay=ay, font=dict(color=INK, size=15),
                bgcolor="rgba(0,0,0,0)", bordercolor=color, borderwidth=1.5, borderpad=3)


# An orthonormal basis Q for C(X), with first column along y_hat. The patch of the plane
# that we draw is a rectangle, in the coordinates of Q, around the points that lie in it.
Q = np.column_stack([unit(y_hat), unit(np.cross(r_hat, y_hat))])
coords = np.vstack([ORIGIN, X.T, y_hat]) @ Q
lo, hi = coords.min(axis=0) - 0.5, coords.max(axis=0) + 0.5
corners = np.array([[lo[0], lo[1]], [hi[0], lo[1]], [hi[0], hi[1]], [lo[0], hi[1]]]) @ Q.T

plane = go.Mesh3d(x=corners[:, 0], y=corners[:, 1], z=corners[:, 2], i=[0, 0], j=[1, 2],
                  k=[2, 3], color=GREEN, opacity=0.18, flatshading=True, hoverinfo="skip",
                  name="𝒞(<b>X</b>), the column space", showlegend=True, legendgroup="plane")

# Faint grid lines in the plane, as a depth cue while rotating.
grid = []
for s in np.linspace(0, 1, 7):
    for a, b in [((lo[0] + s * (hi[0] - lo[0]), lo[1]), (lo[0] + s * (hi[0] - lo[0]), hi[1])),
                 ((lo[0], lo[1] + s * (hi[1] - lo[1])), (hi[0], lo[1] + s * (hi[1] - lo[1])))]:
        grid += [np.array(a) @ Q.T, np.array(b) @ Q.T, np.full(3, np.nan)]
grid = np.array(grid)
plane_grid = go.Scatter3d(x=grid[:, 0], y=grid[:, 1], z=grid[:, 2], mode="lines",
                          line=dict(color=GREEN, width=1.5), opacity=0.45, hoverinfo="skip",
                          showlegend=False, legendgroup="plane")

# The right-angle marker at the foot of the residual: one leg along r_hat, one back along y_hat.
s = 0.3
a, b = s * unit(r_hat), -s * unit(y_hat)
square = np.array([y_hat + a, y_hat + a + b, y_hat + b])
right_angle = go.Scatter3d(x=square[:, 0], y=square[:, 1], z=square[:, 2], mode="lines",
                           line=dict(color=INK, width=4), hoverinfo="skip", showlegend=False,
                           legendgroup="<b>r̂</b> = <b>y</b> − <b>ŷ</b>, the residual")

fig = go.Figure([plane, plane_grid])
fig.add_traces(arrow(ORIGIN, X[:, 0], GREEN, "columns <b>1</b> and <b>x</b> of <b>X</b>",
                     f"<b>1</b> = {vec(X[:, 0])}", width=5, head=0.2, group="cols"))
fig.add_traces(arrow(ORIGIN, X[:, 1], GREEN, "columns <b>1</b> and <b>x</b> of <b>X</b>",
                     f"<b>x</b> = {vec(X[:, 1])}", width=5, head=0.2, group="cols",
                     showlegend=False))
fig.add_traces(arrow(ORIGIN, y, INK, "<b>y</b>, the response", f"<b>y</b> = {vec(y)}"))
fig.add_traces(arrow(ORIGIN, y_hat, BLUE, "<b>ŷ</b> = <b>Py</b>, the fitted values",
                     f"<b>ŷ</b> = <b>X</b>θ̂ = {vec(y_hat)}"
                     f"<br>θ̂ = ({theta_hat[0]:.3f}, {theta_hat[1]:.3f})"))
fig.add_traces(arrow(y_hat, y, ORANGE, "<b>r̂</b> = <b>y</b> − <b>ŷ</b>, the residual",
                     f"<b>r̂</b> = {vec(r_hat)}<br>‖<b>r̂</b>‖ = {np.linalg.norm(r_hat):.3f}"))
fig.add_trace(right_angle)
fig.add_trace(go.Scatter3d(x=[0], y=[0], z=[0], mode="markers", marker=dict(color=INK, size=4),
                           hovertemplate="origin<extra></extra>", showlegend=False))

# Look at the plane obliquely, from the side of it that y is on.
eye = 1.6 * unit(0.4 * Q[:, 0] - 0.6 * Q[:, 1] + unit(r_hat))

# Ranges of equal length on all three axes, so that lengths and angles are not distorted.
pts = np.vstack([ORIGIN, corners, X.T, y])
mid, half = (pts.min(axis=0) + pts.max(axis=0)) / 2, np.ptp(pts, axis=0).max() / 2 + 0.2
axis = dict(showbackground=False, gridcolor=GRID, zerolinecolor=INK, linecolor=GRID,
            tickfont=dict(color=INK))

fig.update_layout(
    template="simple_white",
    paper_bgcolor="rgba(0,0,0,0)",
    font=dict(color=INK, size=13),
    height=640,
    margin=dict(l=0, r=0, t=10, b=0),
    legend=dict(orientation="h", y=0, x=0, yanchor="top"),
    scene=dict(
        xaxis=dict(axis, title="entry 1", range=[mid[0] - half, mid[0] + half]),
        yaxis=dict(axis, title="entry 2", range=[mid[1] - half, mid[1] + half]),
        zaxis=dict(axis, title="entry 3", range=[mid[2] - half, mid[2] + half]),
        aspectmode="manual",
        aspectratio=dict(x=1.25, y=1.25, z=1.25),
        camera=dict(eye=dict(x=eye[0], y=eye[1], z=eye[2]), projection=dict(type="orthographic")),
        annotations=[
            label(y, "<b>y</b>", INK, -35, -25),
            label(0.4 * y_hat, "<b>ŷ</b> = <b>Py</b>", BLUE, -80, 0),
            label(0.5 * (y + y_hat), "<b>r̂</b> = <b>y</b> − <b>ŷ</b>", ORANGE, 0, -60),
            label(X[:, 0], "<b>1</b>", GREEN, 45, 10),
            label(X[:, 1], "<b>x</b>", GREEN, -30, -25),
            label(corners[0] + 0.15 * (corners.mean(axis=0) - corners[0]), "𝒞(<b>X</b>)", GREEN, 0, 0),
            label(ORIGIN, "<b>0</b>", INK, -25, 25),
        ],
    ),
)
fig.show()

Figure 1: Least squares in $\Rbb^3$. The response $\bfy$ (gray) is projected onto the plane $\Ccal(\bfX)$ (green) spanned by the columns $\mathbf{1}$ and $\bfx$. The projection is the vector of fitted values $\widehat{\bfy}$ (blue), and the residual $\widehat{\bfr}=\bfy-\widehat{\bfy}$ (orange) meets the plane at a right angle.

Things to look for while rotating:

$\widehat{\bfy}$ lies in the plane, and is a combination of the two green columns with weights $\widehat{\theta}_1$ and $\widehat{\theta}_2$ (hover over $\widehat{\bfy}$ to see them).
$\widehat{\bfr}$ is perpendicular to both $\widehat{\bfy}$ and the plane. Rotate until you look straight down the edge of the plane and notice that $\widehat{\bfr}$ stands straight up from it.
Any other point $\bfX\bftheta$ of the plane is farther from the tip of $\bfy$, since by the Pythagorean argument in the proof of the normal equations, $\|\bfy-\bfX\bftheta\|_2^2=\|\widehat{\bfr}\|_2^2+\|\bfX(\widehat{\bftheta}-\bftheta)\|_2^2$.

--- title: "The Geometry of Least Squares" subtitle: "Math 123: Mathematical Aspects of Data Analysis — Fall 2026" date: 2026-09-29 number-sections: false jupyter: python3 format: html: toc: true toc-depth: 2 code-line-numbers: true code-tools: true out-width: "100%" column: page execute: echo: true draft: false --- [Least squares is orthogonal projection](../../lectures/lecture-07.qmd#thm-least-squares-projection). That is, the fitted values $\widehat{\bfy}$ are the orthogonal projection of $\bfy$ onto the column space $\Ccal(\bfX)$, and the residual $\widehat{\bfr}=\bfy-\widehat{\bfy}$ is perpendicular to $\Ccal(\bfX)$. ## A tiny regression Consider the setting of fitting a line $y=\theta_1+\theta_2x$ to the three points $(1,1)$, $(2,0)$, $(3,3)$. The design matrix is $\bfX=[\mathbf{1}\ \bfx]\in\Rbb^{3\times2}$, so its columns are the all-ones vector $\mathbf{1}$ and the feature vector $\bfx=(1,2,3)$. ```{python} import numpy as np import plotly.graph_objects as go x = np.array([1.0, 2.0, 3.0]) X = np.column_stack([np.ones(3), x]) y = np.array([1.0, 0.0, 3.0]) theta_hat, *_ = np.linalg.lstsq(X, y, rcond=None) y_hat = X @ theta_hat r_hat = y - y_hat print("theta_hat =", np.round(theta_hat, 4)) print("y_hat =", np.round(y_hat, 4)) print("r_hat =", np.round(r_hat, 4)) print("X^T r_hat =", X.T @ r_hat) ``` The last line is the [normal equations](../../lectures/lecture-07.qmd#thm-normal-equations) $\bfX^\top\widehat{\bfr}=\mathbf{0}$, up to rounding: the residual is orthogonal to both columns of $\bfX$, and hence to the whole plane $\Ccal(\bfX)$. ## The picture In @fig-projection, each axis is one *observation*, not one feature. We are looking at $\Rbb^n$, the space where $\bfy$ lives. Drag to rotate, scroll to zoom, and click legend entries to hide pieces. ```{python} #| label: fig-projection #| fig-cap: "Least squares in $\\Rbb^3$. The response $\\bfy$ (gray) is projected onto the plane $\\Ccal(\\bfX)$ (green) spanned by the columns $\\mathbf{1}$ and $\\bfx$. The projection is the vector of fitted values $\\widehat{\\bfy}$ (blue), and the residual $\\widehat{\\bfr}=\\bfy-\\widehat{\\bfy}$ (orange) meets the plane at a right angle." #| fig-alt: "An interactive three-dimensional plot. A translucent green plane through the origin contains two green arrows, the columns 1 and x of the design matrix. A gray arrow y points from the origin off the plane. A blue arrow y-hat from the origin lies in the plane, directly below the tip of y. An orange arrow r-hat runs from the tip of y-hat to the tip of y, perpendicular to the plane, with a small right-angle marker where it meets y-hat." #| code-fold: true #| code-summary: "Plot the projection geometry" # The same muted palette as the other examples, which reads correctly against # both the light and the dark version of this site. INK = "#7c7c78" GRID = "rgba(128, 128, 128, 0.28)" BLUE, ORANGE, GREEN = "#2a78d6", "#eb6834", "#1baf7a" ORIGIN = np.zeros(3) def unit(v): return v / np.linalg.norm(v) def vec(v): return "(" + ", ".join(f"{c:.3f}" for c in v) + ")" def arrow(start, end, color, name, hover, width=7, head=0.25, showlegend=True, group=None): """A vector from start to end: a line capped by a cone, toggled together in the legend.""" group = group or name d = unit(end - start) shaft = end - 0.5 * head * d # tuck the end of the line inside the cone return [ go.Scatter3d(x=[start[0], shaft[0]], y=[start[1], shaft[1]], z=[start[2], shaft[2]], mode="lines", name=name, legendgroup=group, showlegend=showlegend, line=dict(color=color, width=width), hovertemplate=hover + "<extra></extra>"), go.Cone(x=[end[0]], y=[end[1]], z=[end[2]], u=[head * d[0]], v=[head * d[1]], w=[head * d[2]], anchor="tip", sizemode="raw", sizeref=1, colorscale=[[0, color], [1, color]], showscale=False, lighting=dict(ambient=0.9, diffuse=0.1, specular=0), hoverinfo="skip", legendgroup=group, showlegend=False), ] def label(point, text, color, ax, ay): """A text label pinned to a 3D point, with a thin leader line and a border in the piece's color.""" return dict(x=point[0], y=point[1], z=point[2], text=text, showarrow=True, arrowhead=0, arrowwidth=1, arrowcolor=INK, ax=ax, ay=ay, font=dict(color=INK, size=15), bgcolor="rgba(0,0,0,0)", bordercolor=color, borderwidth=1.5, borderpad=3) # An orthonormal basis Q for C(X), with first column along y_hat. The patch of the plane # that we draw is a rectangle, in the coordinates of Q, around the points that lie in it. Q = np.column_stack([unit(y_hat), unit(np.cross(r_hat, y_hat))]) coords = np.vstack([ORIGIN, X.T, y_hat]) @ Q lo, hi = coords.min(axis=0) - 0.5, coords.max(axis=0) + 0.5 corners = np.array([[lo[0], lo[1]], [hi[0], lo[1]], [hi[0], hi[1]], [lo[0], hi[1]]]) @ Q.T plane = go.Mesh3d(x=corners[:, 0], y=corners[:, 1], z=corners[:, 2], i=[0, 0], j=[1, 2], k=[2, 3], color=GREEN, opacity=0.18, flatshading=True, hoverinfo="skip", name="𝒞(X), the column space", showlegend=True, legendgroup="plane") # Faint grid lines in the plane, as a depth cue while rotating. grid = [] for s in np.linspace(0, 1, 7): for a, b in [((lo[0] + s * (hi[0] - lo[0]), lo[1]), (lo[0] + s * (hi[0] - lo[0]), hi[1])), ((lo[0], lo[1] + s * (hi[1] - lo[1])), (hi[0], lo[1] + s * (hi[1] - lo[1])))]: grid += [np.array(a) @ Q.T, np.array(b) @ Q.T, np.full(3, np.nan)] grid = np.array(grid) plane_grid = go.Scatter3d(x=grid[:, 0], y=grid[:, 1], z=grid[:, 2], mode="lines", line=dict(color=GREEN, width=1.5), opacity=0.45, hoverinfo="skip", showlegend=False, legendgroup="plane") # The right-angle marker at the foot of the residual: one leg along r_hat, one back along y_hat. s = 0.3 a, b = s * unit(r_hat), -s * unit(y_hat) square = np.array([y_hat + a, y_hat + a + b, y_hat + b]) right_angle = go.Scatter3d(x=square[:, 0], y=square[:, 1], z=square[:, 2], mode="lines", line=dict(color=INK, width=4), hoverinfo="skip", showlegend=False, legendgroup="r̂ = y − ŷ, the residual") fig = go.Figure([plane, plane_grid]) fig.add_traces(arrow(ORIGIN, X[:, 0], GREEN, "columns 1 and x of X", f"1 = {vec(X[:, 0])}", width=5, head=0.2, group="cols")) fig.add_traces(arrow(ORIGIN, X[:, 1], GREEN, "columns 1 and x of X", f"x = {vec(X[:, 1])}", width=5, head=0.2, group="cols", showlegend=False)) fig.add_traces(arrow(ORIGIN, y, INK, "y, the response", f"y = {vec(y)}")) fig.add_traces(arrow(ORIGIN, y_hat, BLUE, "ŷ = Py, the fitted values", f"ŷ = Xθ̂ = {vec(y_hat)}" f" θ̂ = ({theta_hat[0]:.3f}, {theta_hat[1]:.3f})")) fig.add_traces(arrow(y_hat, y, ORANGE, "r̂ = y − ŷ, the residual", f"r̂ = {vec(r_hat)} ‖r̂‖ = {np.linalg.norm(r_hat):.3f}")) fig.add_trace(right_angle) fig.add_trace(go.Scatter3d(x=[0], y=[0], z=[0], mode="markers", marker=dict(color=INK, size=4), hovertemplate="origin<extra></extra>", showlegend=False)) # Look at the plane obliquely, from the side of it that y is on. eye = 1.6 * unit(0.4 * Q[:, 0] - 0.6 * Q[:, 1] + unit(r_hat)) # Ranges of equal length on all three axes, so that lengths and angles are not distorted. pts = np.vstack([ORIGIN, corners, X.T, y]) mid, half = (pts.min(axis=0) + pts.max(axis=0)) / 2, np.ptp(pts, axis=0).max() / 2 + 0.2 axis = dict(showbackground=False, gridcolor=GRID, zerolinecolor=INK, linecolor=GRID, tickfont=dict(color=INK)) fig.update_layout( template="simple_white", paper_bgcolor="rgba(0,0,0,0)", font=dict(color=INK, size=13), height=640, margin=dict(l=0, r=0, t=10, b=0), legend=dict(orientation="h", y=0, x=0, yanchor="top"), scene=dict( xaxis=dict(axis, title="entry 1", range=[mid[0] - half, mid[0] + half]), yaxis=dict(axis, title="entry 2", range=[mid[1] - half, mid[1] + half]), zaxis=dict(axis, title="entry 3", range=[mid[2] - half, mid[2] + half]), aspectmode="manual", aspectratio=dict(x=1.25, y=1.25, z=1.25), camera=dict(eye=dict(x=eye[0], y=eye[1], z=eye[2]), projection=dict(type="orthographic")), annotations=[ label(y, "y", INK, -35, -25), label(0.4 * y_hat, "ŷ = Py", BLUE, -80, 0), label(0.5 * (y + y_hat), "r̂ = y − ŷ", ORANGE, 0, -60), label(X[:, 0], "1", GREEN, 45, 10), label(X[:, 1], "x", GREEN, -30, -25), label(corners[0] + 0.15 * (corners.mean(axis=0) - corners[0]), "𝒞(X)", GREEN, 0, 0), label(ORIGIN, "0", INK, -25, 25), ], ), ) fig.show() ``` Things to look for while rotating: - $\widehat{\bfy}$ lies in the plane, and is a combination of the two green columns with weights $\widehat{\theta}_1$ and $\widehat{\theta}_2$ (hover over $\widehat{\bfy}$ to see them). - $\widehat{\bfr}$ is perpendicular to both $\widehat{\bfy}$ *and* the plane. Rotate until you look straight down the edge of the plane and notice that $\widehat{\bfr}$ stands straight up from it. - Any other point $\bfX\bftheta$ of the plane is farther from the tip of $\bfy$, since by the Pythagorean argument in the proof of the [normal equations](../../lectures/lecture-07.qmd#thm-normal-equations), $\|\bfy-\bfX\bftheta\|_2^2=\|\widehat{\bfr}\|_2^2+\|\bfX(\widehat{\bftheta}-\bftheta)\|_2^2$.