SVD: Linear vs. Nonlinear Relationships

Math 123: Mathematical Aspects of Data Analysis — Fall 2026

Published

September 24, 2026

We compare two point clouds in $\Rbb^2$ that are each, up to a little noise, a one-dimensional curve: \[ y=x+\varepsilon \qquad\text{and}\qquad y=x^2+\varepsilon, \] with $x$ uniform on $[-1,1]$ and $\varepsilon$ small Gaussian noise. In both cases $y$ is (almost) a function of $x$, so each cloud really has one degree of freedom. The SVD finds this for the line and misses it for the parabola, because the best it can do is a best-fitting subspace, and a parabola is not close to any line.

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

rng = np.random.default_rng(123)
n, noise = 200, 0.05
x = rng.uniform(-1, 1, n)
z = rng.standard_normal(n)             # one fixed draw of N(0, 1), rescaled below
eps = noise * z

relations = {"linear": lambda t: t, "quadratic": lambda t: t**2}
data = {name: np.column_stack([x, f(x) + eps]) for name, f in relations.items()}


def centered_svd(D):
    """Center the rows of D and return the mean, the centered matrix X, and the SVD of X."""
    mean = D.mean(axis=0)
    X = D - mean
    U, s, Vt = np.linalg.svd(X, full_matrices=False)
    V = Vt.T
    # Singular vectors are only determined up to sign; make v_1 point right and v_2 point up.
    for j, flip in enumerate([V[0, 0] < 0, V[1, 1] < 0]):
        if flip:
            V[:, j], U[:, j] = -V[:, j], -U[:, j]
    return mean, X, U, s, V


fits = {name: centered_svd(D) for name, D in data.items()}
for name, (mean, X, U, s, V) in fits.items():
    R1 = s[0] ** 2 / np.sum(s**2)
    print(f"{name:>9}:  sigma/sqrt(n) = {np.round(s / np.sqrt(n), 4)},  "
          f"R_1 = {R1:.4f},  v_1 = {np.round(V[:, 0], 3)}")

   linear:  sigma/sqrt(n) = [0.8096 0.0387],  R_1 = 0.9977,  v_1 = [0.707 0.707]
quadratic:  sigma/sqrt(n) = [0.5733 0.2955],  R_1 = 0.7901,  v_1 = [ 0.999 -0.042]

Figure style and plotting helpers

# 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 = "#2a78d6", "#eb6834"
TITLES = ["y = x + ε", "y = x² + ε"]


def style(fig, height=460, title=None, showlegend=None):
    """Apply the shared layout to a plot with axes."""
    fig.update_layout(
        template="simple_white",
        paper_bgcolor="rgba(0,0,0,0)",
        plot_bgcolor="rgba(0,0,0,0)",
        font=dict(color=INK, size=13),
        title=dict(text=title, x=0, xanchor="left", font=dict(size=15)),
        margin=dict(l=60, r=30, t=60 if title else 40, b=50),
        height=height,
        showlegend=showlegend,
        hovermode="closest",
    )
    fig.update_xaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID)
    fig.update_yaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID)
    return fig


def equal_axes(fig, lim, row=None, col=None, anchor="x"):
    """Square axes on [-lim, lim], so right angles in the data look like right angles."""
    fig.update_xaxes(range=[-lim, lim], constrain="domain", row=row, col=col)
    fig.update_yaxes(range=[-lim, lim], scaleanchor=anchor, scaleratio=1, row=row, col=col)


def points(P, name, color=INK, size=6, opacity=0.7, showlegend=True):
    """Markers at the rows of P."""
    return go.Scatter(x=P[:, 0], y=P[:, 1], mode="markers", name=name, legendgroup=name,
                      showlegend=showlegend,
                      hovertemplate="(%{x:.3f}, %{y:.3f})<extra>" + name + "</extra>",
                      marker=dict(color=color, size=size, opacity=opacity,
                                  line=dict(color="white", width=0.5)))


def arrow(start, tip, color, name, showlegend=True):
    """An arrow from start to tip."""
    return go.Scatter(
        x=[start[0], tip[0]], y=[start[1], tip[1]], mode="lines+markers", name=name,
        legendgroup=name, showlegend=showlegend, hoverinfo="skip",
        line=dict(color=color, width=3),
        marker=dict(symbol="arrow", angleref="previous", size=[0, 14], color=color),
    )


def segments(P, Q, name, showlegend=True):
    """Thin segments from each row of P to the matching row of Q, drawn as one trace."""
    xs = np.column_stack([P[:, 0], Q[:, 0], np.full(len(P), np.nan)]).ravel()
    ys = np.column_stack([P[:, 1], Q[:, 1], np.full(len(P), np.nan)]).ravel()
    return go.Scatter(x=xs, y=ys, mode="lines", name=name, legendgroup=name,
                      showlegend=showlegend, hoverinfo="skip",
                      line=dict(color=GRID, width=1))

The principal directions

Figure 1 shows both clouds together with their principal directions $\bfv_1,\bfv_2$, drawn from the mean $\bar{\bfx}$. Each arrow has length $2\sigma_j/\sqrt n$, twice the standard deviation of the data along $\bfv_j$, so it shows how far the cloud spreads in that direction.

Plot the data and the principal directions

fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.1, subplot_titles=TITLES)
for col, (name, (mean, X, U, s, V)) in enumerate(fits.items(), start=1):
    first = col == 1
    fig.add_trace(points(data[name], "data", showlegend=first), row=1, col=col)
    fig.add_trace(arrow(mean, mean + 2 * s[0] / np.sqrt(n) * V[:, 0], BLUE,
                        "2(σ<sub>1</sub>/√n) v<sub>1</sub>", showlegend=first), row=1, col=col)
    fig.add_trace(arrow(mean, mean + 2 * s[1] / np.sqrt(n) * V[:, 1], ORANGE,
                        "2(σ<sub>2</sub>/√n) v<sub>2</sub>", showlegend=first), row=1, col=col)
    equal_axes(fig, 1.3, row=1, col=col, anchor="x" if first else "x2")
style(fig, height=500, showlegend=True)
fig.update_layout(legend=dict(orientation="h", y=-0.12, x=0.3))
fig.show()

Figure 1: The two point clouds, each with arrows 2 (sigma_j / sqrt(n)) v_j drawn from the sample mean.

For the line, $\sigma_2$ is tiny and $\bfv_1\approx(1,1)/\sqrt2$ points along the line, so the explained variance ratio $R_1$ is essentially $1$. For the parabola, $\bfv_1\approx(1,0)$ is just the $x$-axis, $\bfv_2\approx(0,1)$ is the $y$-axis, and $R_1$ is only about $0.79$. As far as the SVD is concerned, the parabola is a genuinely two-dimensional cloud.

The rank-1 approximation

By best-fitting subspaces come from the SVD, the rank-1 reconstructions $\widehat{\bfx}_i=\bar{\bfx}+\bfv_1\bfv_1^\top(\bfx_i-\bar{\bfx})$ are the orthogonal projections of the data onto the line through $\bar{\bfx}$ in the direction $\bfv_1$. Figure 2 draws them, with a thin segment from each data point to its reconstruction.

Plot the rank-1 reconstructions

fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.1, subplot_titles=TITLES)
for col, (name, (mean, X, U, s, V)) in enumerate(fits.items(), start=1):
    first = col == 1
    recon = mean + np.outer(X @ V[:, 0], V[:, 0])      # rows are x_bar + v_1 v_1^T (x_i - x_bar)
    fig.add_trace(segments(data[name], recon, "residual", showlegend=first), row=1, col=col)
    fig.add_trace(points(data[name], "data", opacity=0.5, showlegend=first), row=1, col=col)
    fig.add_trace(points(recon, "rank-1 reconstruction", color=BLUE, size=6, opacity=0.9,
                         showlegend=first), row=1, col=col)
    equal_axes(fig, 1.3, row=1, col=col, anchor="x" if first else "x2")
style(fig, height=500, showlegend=True)
fig.update_layout(legend=dict(orientation="h", y=-0.12, x=0.25))
fig.show()

Figure 2: The rank-1 reconstructions (blue) of each point cloud (grey), with a segment from each point to its reconstruction. The line survives; the parabola is flattened onto the horizontal line y = mean of y.

For the parabola, the rank-1 reconstruction throws away $y$ entirely: every point is sent to $(x_i,\bar y)$. This is despite $y$ being determined by $x$ up to noise, so a single number per point really would suffice to describe the data. The SVD just cannot express $y=x^2$ as, fundamentally, this is a linear dimensionality reduction tool.

Hope on the Horizon

The relationship $y=x^2$ is nonlinear in $x$, but it is linear in the features $(x,x^2)$. If we add $x^2$ as a third column, the data $(x_i,x_i^2,y_i)\in\Rbb^3$ lie near the plane $\{(a,b,c):c=b\}$, and the SVD finds it: $\sigma_3/\sqrt n$ drops to the size of the noise, just like $\sigma_2/\sqrt n$ for the line. Choosing good features, by hand or implicitly via kernel methods, is one standard way to let linear tools like the SVD see nonlinear structure.

lifted = np.column_stack([x, x**2, data["quadratic"][:, 1]])
s_lifted = np.linalg.svd(lifted - lifted.mean(axis=0), compute_uv=False)
print("features (x, y):      sigma/sqrt(n) =", np.round(fits["quadratic"][3] / np.sqrt(n), 4))
print("features (x, x^2, y): sigma/sqrt(n) =", np.round(s_lifted / np.sqrt(n), 4))
print(f"noise / sqrt(2) = {noise / np.sqrt(2):.4f}")

features (x, y):      sigma/sqrt(n) = [0.5733 0.2955]
features (x, x^2, y): sigma/sqrt(n) = [0.5739 0.4086 0.038 ]
noise / sqrt(2) = 0.0354

--- title: "SVD: Linear vs. Nonlinear Relationships" subtitle: "Math 123: Mathematical Aspects of Data Analysis — Fall 2026" date: 2026-09-24 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 --- We compare two point clouds in $\Rbb^2$ that are each, up to a little noise, a one-dimensional curve: $$ y=x+\varepsilon \qquad\text{and}\qquad y=x^2+\varepsilon, $$ with $x$ uniform on $[-1,1]$ and $\varepsilon$ small Gaussian noise. In both cases $y$ is (almost) a function of $x$, so each cloud really has one degree of freedom. The SVD finds this for the line and misses it for the parabola, because the best it can do is a [best-fitting subspace](../../lectures/lecture-06.qmd#thm-best-subspace), and a parabola is not close to any line. ```{python} import numpy as np import plotly.graph_objects as go from plotly.subplots import make_subplots rng = np.random.default_rng(123) n, noise = 200, 0.05 x = rng.uniform(-1, 1, n) z = rng.standard_normal(n) # one fixed draw of N(0, 1), rescaled below eps = noise * z relations = {"linear": lambda t: t, "quadratic": lambda t: t**2} data = {name: np.column_stack([x, f(x) + eps]) for name, f in relations.items()} def centered_svd(D): """Center the rows of D and return the mean, the centered matrix X, and the SVD of X.""" mean = D.mean(axis=0) X = D - mean U, s, Vt = np.linalg.svd(X, full_matrices=False) V = Vt.T # Singular vectors are only determined up to sign; make v_1 point right and v_2 point up. for j, flip in enumerate([V[0, 0] < 0, V[1, 1] < 0]): if flip: V[:, j], U[:, j] = -V[:, j], -U[:, j] return mean, X, U, s, V fits = {name: centered_svd(D) for name, D in data.items()} for name, (mean, X, U, s, V) in fits.items(): R1 = s[0] ** 2 / np.sum(s**2) print(f"{name:>9}: sigma/sqrt(n) = {np.round(s / np.sqrt(n), 4)}, " f"R_1 = {R1:.4f}, v_1 = {np.round(V[:, 0], 3)}") ``` ```{python} #| code-fold: true #| code-summary: "Figure style and plotting helpers" # 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 = "#2a78d6", "#eb6834" TITLES = ["y = x + ε", "y = x² + ε"] def style(fig, height=460, title=None, showlegend=None): """Apply the shared layout to a plot with axes.""" fig.update_layout( template="simple_white", paper_bgcolor="rgba(0,0,0,0)", plot_bgcolor="rgba(0,0,0,0)", font=dict(color=INK, size=13), title=dict(text=title, x=0, xanchor="left", font=dict(size=15)), margin=dict(l=60, r=30, t=60 if title else 40, b=50), height=height, showlegend=showlegend, hovermode="closest", ) fig.update_xaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID) fig.update_yaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID) return fig def equal_axes(fig, lim, row=None, col=None, anchor="x"): """Square axes on [-lim, lim], so right angles in the data look like right angles.""" fig.update_xaxes(range=[-lim, lim], constrain="domain", row=row, col=col) fig.update_yaxes(range=[-lim, lim], scaleanchor=anchor, scaleratio=1, row=row, col=col) def points(P, name, color=INK, size=6, opacity=0.7, showlegend=True): """Markers at the rows of P.""" return go.Scatter(x=P[:, 0], y=P[:, 1], mode="markers", name=name, legendgroup=name, showlegend=showlegend, hovertemplate="(%{x:.3f}, %{y:.3f})<extra>" + name + "</extra>", marker=dict(color=color, size=size, opacity=opacity, line=dict(color="white", width=0.5))) def arrow(start, tip, color, name, showlegend=True): """An arrow from start to tip.""" return go.Scatter( x=[start[0], tip[0]], y=[start[1], tip[1]], mode="lines+markers", name=name, legendgroup=name, showlegend=showlegend, hoverinfo="skip", line=dict(color=color, width=3), marker=dict(symbol="arrow", angleref="previous", size=[0, 14], color=color), ) def segments(P, Q, name, showlegend=True): """Thin segments from each row of P to the matching row of Q, drawn as one trace.""" xs = np.column_stack([P[:, 0], Q[:, 0], np.full(len(P), np.nan)]).ravel() ys = np.column_stack([P[:, 1], Q[:, 1], np.full(len(P), np.nan)]).ravel() return go.Scatter(x=xs, y=ys, mode="lines", name=name, legendgroup=name, showlegend=showlegend, hoverinfo="skip", line=dict(color=GRID, width=1)) ``` ## The principal directions {#sec-directions} @fig-data shows both clouds together with their [principal directions](../../lectures/lecture-06.qmd#def-pca) $\bfv_1,\bfv_2$, drawn from the mean $\bar{\bfx}$. Each arrow has length $2\sigma_j/\sqrt n$, twice the standard deviation of the data along $\bfv_j$, so it shows how far the cloud spreads in that direction. ```{python} #| label: fig-data #| fig-cap: "The two point clouds, each with arrows 2 (sigma_j / sqrt(n)) v_j drawn from the sample mean." #| fig-alt: "Two square scatter plots of 200 grey points. On the left the points lie along the diagonal line y = x, with a long blue arrow along the diagonal and a barely visible orange arrow perpendicular to it. On the right the points lie along an upward parabola, with a long blue arrow pointing horizontally and a shorter but clearly visible orange arrow pointing vertically." #| code-fold: true #| code-summary: "Plot the data and the principal directions" fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.1, subplot_titles=TITLES) for col, (name, (mean, X, U, s, V)) in enumerate(fits.items(), start=1): first = col == 1 fig.add_trace(points(data[name], "data", showlegend=first), row=1, col=col) fig.add_trace(arrow(mean, mean + 2 * s[0] / np.sqrt(n) * V[:, 0], BLUE, "2(σ1/√n) v1", showlegend=first), row=1, col=col) fig.add_trace(arrow(mean, mean + 2 * s[1] / np.sqrt(n) * V[:, 1], ORANGE, "2(σ2/√n) v2", showlegend=first), row=1, col=col) equal_axes(fig, 1.3, row=1, col=col, anchor="x" if first else "x2") style(fig, height=500, showlegend=True) fig.update_layout(legend=dict(orientation="h", y=-0.12, x=0.3)) fig.show() ``` For the line, $\sigma_2$ is tiny and $\bfv_1\approx(1,1)/\sqrt2$ points along the line, so the [explained variance ratio](../../lectures/lecture-06.qmd#def-explained-variance) $R_1$ is essentially $1$. For the parabola, $\bfv_1\approx(1,0)$ is just the $x$-axis, $\bfv_2\approx(0,1)$ is the $y$-axis, and $R_1$ is only about $0.79$. As far as the SVD is concerned, the parabola is a genuinely two-dimensional cloud. ## The rank-1 approximation {#sec-rank-one} By [best-fitting subspaces come from the SVD](../../lectures/lecture-06.qmd#thm-best-subspace), the rank-1 reconstructions $\widehat{\bfx}_i=\bar{\bfx}+\bfv_1\bfv_1^\top(\bfx_i-\bar{\bfx})$ are the orthogonal projections of the data onto the line through $\bar{\bfx}$ in the direction $\bfv_1$. @fig-rank-one draws them, with a thin segment from each data point to its reconstruction. ```{python} #| label: fig-rank-one #| fig-cap: "The rank-1 reconstructions (blue) of each point cloud (grey), with a segment from each point to its reconstruction. The line survives; the parabola is flattened onto the horizontal line y = mean of y." #| fig-alt: "Two square panels. On the left the blue reconstructed points lie on the diagonal line almost exactly on top of the grey data, with invisible connecting segments. On the right the grey points form an upward parabola while the blue reconstructed points all lie on a horizontal line at height about one third, joined to the grey points by long vertical segments." #| code-fold: true #| code-summary: "Plot the rank-1 reconstructions" fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.1, subplot_titles=TITLES) for col, (name, (mean, X, U, s, V)) in enumerate(fits.items(), start=1): first = col == 1 recon = mean + np.outer(X @ V[:, 0], V[:, 0]) # rows are x_bar + v_1 v_1^T (x_i - x_bar) fig.add_trace(segments(data[name], recon, "residual", showlegend=first), row=1, col=col) fig.add_trace(points(data[name], "data", opacity=0.5, showlegend=first), row=1, col=col) fig.add_trace(points(recon, "rank-1 reconstruction", color=BLUE, size=6, opacity=0.9, showlegend=first), row=1, col=col) equal_axes(fig, 1.3, row=1, col=col, anchor="x" if first else "x2") style(fig, height=500, showlegend=True) fig.update_layout(legend=dict(orientation="h", y=-0.12, x=0.25)) fig.show() ``` For the parabola, the rank-1 reconstruction throws away $y$ entirely: every point is sent to $(x_i,\bar y)$. This is despite $y$ being determined by $x$ up to noise, so a single number per point really would suffice to describe the data. The SVD just cannot express $y=x^2$ as, fundamentally, this is a linear dimensionality reduction tool. ::: {.callout-tip} ## Hope on the Horizon The relationship $y=x^2$ is nonlinear in $x$, but it is *linear* in the features $(x,x^2)$. If we add $x^2$ as a third column, the data $(x_i,x_i^2,y_i)\in\Rbb^3$ lie near the plane $\{(a,b,c):c=b\}$, and the SVD finds it: $\sigma_3/\sqrt n$ drops to the size of the noise, just like $\sigma_2/\sqrt n$ for the line. Choosing good features, by hand or implicitly via kernel methods, is one standard way to let linear tools like the SVD see nonlinear structure. ```{python} lifted = np.column_stack([x, x**2, data["quadratic"][:, 1]]) s_lifted = np.linalg.svd(lifted - lifted.mean(axis=0), compute_uv=False) print("features (x, y): sigma/sqrt(n) =", np.round(fits["quadratic"][3] / np.sqrt(n), 4)) print("features (x, x^2, y): sigma/sqrt(n) =", np.round(s_lifted / np.sqrt(n), 4)) print(f"noise / sqrt(2) = {noise / np.sqrt(2):.4f}") ``` :::