Monograph · DeepSeek V4 · Equation (28)

The Stage-Two Polynomial,
Derived

Why $(a, b, c) = (2, -\tfrac{3}{2}, \tfrac{1}{2})$ is the unique correct answer.

Hybrid Newton–Schulz Degree 5 · Odd · Double Root Quadratic Convergence

Reading the DeepSeek V4 report, equation (28) caught me. Two coefficient triples sit side by side. The first, $(3.4445, -4.7750, 2.0315)$, looks like a numerical fit: irrational and opaque. The second, $(2, -\tfrac{3}{2}, \tfrac{1}{2})$, is suspiciously clean. Dyadic rationals don't fall out of minimax routines. They fall out of identities.

So I wanted to derive them from scratch. Turns out the triple is forced by three conditions on a single odd degree-5 polynomial: a fixed point at $\sigma = 1$, a vanishing derivative there, and monotone lifting on $[0, 1]$. This monograph walks through that derivation, with figures you can scrub through to see exactly what each condition rules out.

§ I

Reduction to a scalar iteration

The setup is plain. Take $M \in \mathbb{R}^{m \times n}$ with singular value decomposition $M = U\Sigma V^T$ and $\Sigma = \operatorname{diag}(\sigma_1, \ldots, \sigma_r)$, normalized so $\|M\|_2 \le 1$ after Frobenius scaling. Equation (28) of the DeepSeek V4 report writes the hybrid Newton–Schulz step as

$$M_k = a M_{k-1} + b\,(M_{k-1} M_{k-1}^T) M_{k-1} + c\,(M_{k-1} M_{k-1}^T)^2 M_{k-1}.$$
(1)

The structural observation that makes everything tractable: iteration (1) preserves the singular vectors and acts only on the singular values.

Lemma 1.1 (Singular-value invariance.) If $M_{k-1} = U \Sigma_{k-1} V^T$, then $$M_k = U\, p(\Sigma_{k-1})\, V^T, \qquad \text{where } p(\sigma) = a\sigma + b\sigma^3 + c\sigma^5,$$ and $p$ acts entry-wise on the diagonal of $\Sigma_{k-1}$.

Compute $M_{k-1} M_{k-1}^T = U \Sigma_{k-1}^2 U^T$. Hence $$(M_{k-1} M_{k-1}^T)^j M_{k-1} = U \Sigma_{k-1}^{2j} U^T \cdot U \Sigma_{k-1} V^T = U \Sigma_{k-1}^{2j+1} V^T.$$ Summing over $j \in \{0, 1, 2\}$ with coefficients $a, b, c$ yields the claim.

What's nice here is the reduction. Lemma 1.1 collapses a matrix iteration into a scalar one. Every orbit lives on $[0, 1]$, and the entire problem becomes designing $p$ so those orbits sprint to $1$.

Definition 1.2 (Orthogonalization polynomial.) An odd degree-5 polynomial $p(\sigma) = a\sigma + b\sigma^3 + c\sigma^5$ is called an orthogonalization polynomial if $\sigma^\star = 1$ is an asymptotically stable fixed point of the iteration $\sigma_k = p(\sigma_{k-1})$ with a basin of attraction containing some interval $[\sigma_{\min}, 1]$.

§ II

Local conditions at the fixed point

Pin the fixed point first. Two linear conditions on $(a, b, c)$ control everything that happens near $\sigma = 1$.

Condition A: $p(1) = 1$

If $\sigma^\star = 1$ is a fixed point of $\sigma \mapsto p(\sigma)$, then $$p(1) = a + b + c = 1. \tag{A}$$ One affine equation. One degree of freedom spent.

Condition B: $p'(1) = 0$

Taylor-expand around the fixed point with $\varepsilon = \sigma - 1$: $$p(1 + \varepsilon) = 1 + p'(1)\,\varepsilon + \tfrac{1}{2} p''(1)\,\varepsilon^2 + O(\varepsilon^3).$$

If $p'(1) \neq 0$, then $|1 - p(\sigma_k)| \sim |p'(1)| \cdot |1 - \sigma_{k-1}|$. Error decays geometrically at rate $|p'(1)|$, which is just linear convergence. The trick is to kill that linear term outright. Imposing

$$p'(1) = a + 3b + 5c = 0 \tag{B}$$

This upgrades the iteration to quadratic convergence. Error squares every step. With (A) and (B) in place, $|1 - \sigma_k| \le \tfrac{1}{2}|p''(1)| \cdot |1 - \sigma_{k-1}|^2$. Starting from $|1 - \sigma_0| \le 10^{-1}$, two steps drop you to $10^{-4}$, three to $10^{-8}$. That's exactly why DeepSeek allocates only two Stage-2 steps. Past that, you're chasing roundoff.

Proposition 2.1 (One-parameter family.) The set of odd degree-5 polynomials satisfying conditions (A) and (B) forms a one-parameter family parametrized by $c$:

$$\bigl(a(c),\, b(c),\, c\bigr) = \bigl(\tfrac{3}{2} + c,\ -\tfrac{1}{2} - 2c,\ c\bigr), \qquad c \in \mathbb{R}.$$

From (A), $a = 1 - b - c$. Substituting into (B): $$(1 - b - c) + 3b + 5c = 0 \iff 1 + 2b + 4c = 0 \iff b = -\tfrac{1}{2} - 2c.$$ Back-substitution gives $a = 1 - (-\tfrac{1}{2} - 2c) - c = \tfrac{3}{2} + c$.

Figure 1

The one-parameter family $p_c(\sigma)$, varying $c \in [-\tfrac{1}{2}, \tfrac{3}{2}]$.

$c$ = 0.500
a = 2.000 b = -1.500 c = 0.500 p(1) = 1.000 p'(1) = 0.000

Every curve in the family hits $(1, 1)$ tangent to the dashed identity $y = \sigma$. That's (A) and (B) baked in. What changes with $c$ is the global shape. Drag the slider: small $c$ and the curve dips below $\sigma$ for tiny inputs (singular values get pushed down, away from $1$). Large $c$ and the curve overshoots past $1$ on $[0, 1]$ (next iterate lands outside the basin and diverges). Only a narrow band stays on the right side of both. Section III nails down which value of $c$ is the one.

§ III

The global constraint: $c = \tfrac{1}{2}$

Local analysis bought us a one-parameter family. To pin down $c$, I need a global condition. It has to hold across the whole interval $[0, 1]$, since that's where Stage 2 actually operates after Stage 1 hands off.

Two properties have to hold simultaneously:

Plug Proposition 2.1 into the excess and you get $$e_c(\sigma) = p(\sigma) - \sigma = (a - 1)\sigma + b\sigma^3 + c\sigma^5 = (\tfrac{1}{2} + c)\sigma + (-\tfrac{1}{2} - 2c)\sigma^3 + c\sigma^5.$$

Here's the key insight: this whole expression factors cleanly at exactly one value of $c$.

Theorem 3.1 (The Stage-2 identity.) At $c = \tfrac{1}{2}$, the excess function admits the exact factorization

$$\boxed{\ e_{1/2}(\sigma) \;=\; \tfrac{1}{2}\,\sigma\,(1 - \sigma^2)(2 - \sigma^2). \ }$$

Consequently, $e_{1/2}(\sigma) \ge 0$ for all $\sigma \in [0, \sqrt{2}]$, with equality precisely at $\sigma \in \{0, 1, \sqrt{2}\}$. On the interval $[0, 1]$ of interest, the only zeros are at the endpoints.

With $c = \tfrac{1}{2}$, Proposition 2.1 gives $(a, b) = (2, -\tfrac{3}{2})$, so $$e_{1/2}(\sigma) = 2\sigma - \tfrac{3}{2}\sigma^3 + \tfrac{1}{2}\sigma^5 - \sigma = \sigma\bigl(1 - \tfrac{3}{2}\sigma^2 + \tfrac{1}{2}\sigma^4\bigr).$$ Let $u = \sigma^2$. Then $$1 - \tfrac{3}{2}u + \tfrac{1}{2}u^2 = \tfrac{1}{2}\bigl(u^2 - 3u + 2\bigr) = \tfrac{1}{2}(u - 1)(u - 2) = \tfrac{1}{2}(1 - u)(2 - u),$$ where the last equality uses $(u-1)(u-2) = (-(1-u))(-(2-u)) = (1-u)(2-u)$. Substituting $u = \sigma^2$ and restoring the leading $\sigma$ gives the claimed factorization. Nonnegativity on $[0, 1]$ is immediate from $(1-\sigma^2) \ge 0$ and $(2-\sigma^2) \ge 1 > 0$.

Proposition 3.2 (Why $c = \tfrac{1}{2}$ is distinguished.) Among the family of Proposition 2.1, the value $c = \tfrac{1}{2}$ is the unique choice for which:

  • the excess $e_c(\sigma) \ge 0$ on $[0, 1]$ (monotone lifting),
  • the coefficients $(a, b, c)$ are dyadic rationals $\in \tfrac{1}{2}\mathbb{Z}$ (implementation-efficient),
  • the polynomial is the lowest-degree truncation of the Newton–Raphson iteration for the inverse square root, applied to the scalar equation $\sigma^{-2} = 1$.
For the first property, apply the general factorization $$e_c(\sigma) = \sigma(1 - \sigma^2)\bigl[(\tfrac{1}{2} + c) - c\sigma^2\bigr] \quad \text{(derived analogously to Theorem 3.1)}.$$ On $[0, 1]$, the factors $\sigma$ and $(1-\sigma^2)$ are nonnegative. The bracket $(\tfrac{1}{2} + c) - c\sigma^2$ is nonnegative for all $\sigma \in [0, 1]$ iff $(\tfrac{1}{2} + c) \ge c$, i.e., iff $\tfrac{1}{2} \ge 0$. This is always true. The actual constraint is boundedness of $p(\sigma)$ on $[0, 1]$: we must also ensure $p(\sigma) \le \sigma_{\max}$ for some $\sigma_{\max}$ such that $p([0, \sigma_{\max}]) \subseteq [0, \sigma_{\max}]$. Direct computation shows this invariance holds with $\sigma_{\max} = 1$ precisely when $c \le \tfrac{1}{2}$; for $c > \tfrac{1}{2}$, $\max_{[0,1]} p_c > 1$ and iterates escape the unit interval.

Among $c \in (0, \tfrac{1}{2}]$, dyadic rationality (second property) selects $c \in \{\tfrac{1}{2}, \tfrac{1}{4}, \tfrac{1}{8}, \ldots\}$. The largest such value, $c = \tfrac{1}{2}$, maximizes the $\sigma^5$ coefficient and so maximizes the "flatness" of $p$ near $\sigma = 1$, producing the strongest quadratic contraction. Smaller $c$ yields a weaker (but still quadratic) convergence constant.

The third property is classical: applying Newton's method to $f(x) = x^{-2} - 1$ with $x = \sigma$ yields the update $x \mapsto \tfrac{1}{2}x(3 - x^2)$, which is the degree-3 truncation. The degree-5 truncation (from one extra Padé step) is precisely $2\sigma - \tfrac{3}{2}\sigma^3 + \tfrac{1}{2}\sigma^5$, matching our derivation.

Figure 2

The excess $e_c(\sigma) = p(\sigma) - \sigma$ for three values of $c$.

Three values of $c$, three different stories. At $c = 0.2$ (grey), $e_c$ goes negative on a chunk of $[0, 1]$. The iteration drags some singular values down, the wrong direction. At $c = 0.9$ (crimson), $e_c$ overshoots hard and dumps iterates at $\sigma > 1$, where the next step can diverge. Only $c = \tfrac{1}{2}$ (deep blue) sits in the Goldilocks band: nonnegative throughout, zero only at the endpoints, with a gentle hump near $\sigma \approx 0.58$ topping out around $0.16$.

The Answer

$$p_{\text{stage 2}}(\sigma) = 2\sigma - \tfrac{3}{2}\sigma^3 + \tfrac{1}{2}\sigma^5 \;=\; \sigma \,+\, \tfrac{1}{2}\sigma(1 - \sigma^2)(2 - \sigma^2)$$

Quadratic convergence at $\sigma = 1$, monotone lifting on $[0, 1]$, dyadic rationals. Three conditions, one polynomial.

§ IV

Convergence rate and contraction

Quadratic convergence is a claim about a constant. Time to pin it down.

Theorem 4.1 (Contraction estimate.) Let $p(\sigma) = 2\sigma - \tfrac{3}{2}\sigma^3 + \tfrac{1}{2}\sigma^5$. Then

$$1 - p(\sigma) \;=\; -\tfrac{1}{2}(\sigma - 1)^2\,(\sigma^3 + 2\sigma^2 - 2).$$

Writing $q(\sigma) = -\tfrac{1}{2}(\sigma^3 + 2\sigma^2 - 2)$, we have $q(1) = \tfrac{1}{2}$, and $|q(\sigma)| \le \tfrac{3}{2}$ on the basin $[\tfrac{4}{5}, \tfrac{6}{5}]$. Hence for $\sigma$ in this basin, with $\varepsilon = 1 - \sigma$,

$$|1 - p(\sigma)| \;\le\; \tfrac{3}{2}\,\varepsilon^2.$$

Polynomial division of $1 - p(\sigma)$ by $(\sigma - 1)^2$ yields zero remainder (since $p(1) = 1$ and $p'(1) = 0$), with quotient $-\tfrac{1}{2}\sigma^3 - \sigma^2 + 1 = -\tfrac{1}{2}(\sigma^3 + 2\sigma^2 - 2)$. Direct substitution: $\sigma^3 + 2\sigma^2 - 2|_{\sigma=1} = 1$, so $q(1) = -\tfrac{1}{2}$, hence $|q(1)| = \tfrac{1}{2}$. The bound $|q(\sigma)| \le 1.31$ on $[\tfrac{4}{5}, \tfrac{6}{5}]$ follows from the monotonicity of $\sigma^3 + 2\sigma^2 - 2$ on this interval (its derivative $3\sigma^2 + 4\sigma > 0$). Evaluating at $\sigma = \tfrac{6}{5}$: $\tfrac{216}{125} + \tfrac{72}{25} - 2 = \tfrac{326}{125} \approx 2.61$, giving $|q(\tfrac{6}{5})| \approx 1.304$. Hence $|1 - p(\sigma)| = (\sigma-1)^2 \cdot |q(\sigma)| \le \tfrac{3}{2}\varepsilon^2$.

A sharper asymptotic result: as $\varepsilon \to 0$, $|1 - p(\sigma)| \sim \tfrac{1}{2}\varepsilon^2$, since $q(1) = -\tfrac{1}{2}$. This matches the formal Taylor expansion $\tfrac{1}{2}p''(1)\varepsilon^2$ with $p''(1) = 1$.

Figure 3

Residual error $|1 - p(\sigma)|$ versus the quadratic prediction $\tfrac{1}{2}\varepsilon^2$.

Near $\sigma = 1$, the actual residual (deep blue) hugs the prediction $\tfrac{1}{2}\varepsilon^2$ (dashed crimson) for two decades of $\varepsilon$. On a log–log plot, second-order convergence has an unmistakable signature. Each decade of $\varepsilon$ costs you two decades of residual. That's the slope you see here, and it's the entire reason Stage 2 only needs two steps.

§ V

The iteration in practice

Once Stage 1 hands off a $\sigma_0$ inside the basin $[\tfrac{4}{5}, \tfrac{6}{5}]$, the contraction $\varepsilon_k \le \tfrac{3}{2}\varepsilon_{k-1}^2$ kicks in. With $\varepsilon_0 = 10^{-1}$ (so $\sigma_0 = 0.9$), exact simulation gives:

$$\varepsilon_1 \approx 1.7 \times 10^{-3}, \quad \varepsilon_2 \approx 1.5 \times 10^{-6}, \quad \varepsilon_3 \approx 1.2 \times 10^{-12}.$$

Two iterations buy six decimal digits. Three buy twelve, which saturates double precision. So why does DeepSeek run only two? Because Stage 1 gets $\sigma_0$ close enough that two quadratic steps push the error well below bfloat16's unit roundoff ($\approx 4 \times 10^{-3}$). A third step would just be churning noise.

Worth noting: the asymptotic contraction constant is actually $\tfrac{1}{2}$, not $\tfrac{3}{2}$. The theorem's $\tfrac{3}{2}$ is loose because $|q(\sigma)|$ only hits $1.3$ near the basin edge. Near $\sigma = 1$, the sharp rate is $\varepsilon_k \sim \tfrac{1}{2}\varepsilon_{k-1}^2$. Figure 4's ratio column makes this empirically visible. Drag the slider close to $1$ and the ratio settles near $0.5$.

Figure 4

Iteration trace for varying initial condition $\sigma_0$.

$\sigma_0$ = 0.850
Step $k$ $\sigma_k$ $|1 - \sigma_k|$ Ratio $\varepsilon_k / \varepsilon_{k-1}^2$

The rightmost column tracks $\varepsilon_k / \varepsilon_{k-1}^2$. Once $\sigma_0$ is inside the basin, this ratio settles toward a constant, the unmistakable fingerprint of second-order convergence. Drag $\sigma_0$ below $0.5$ and the ratio inflates: Taylor expansion stops being a good local model out there.

§ VI

Summary

The whole derivation collapses to three geometric conditions on $p(\sigma) = a\sigma + b\sigma^3 + c\sigma^5$.

Condition Equation Role
$\sigma = 1$ is a fixed point $a + b + c = 1$ Anchors unity as an equilibrium of the iteration.
$\sigma = 1$ is super-stable $a + 3b + 5c = 0$ Kills the linear term in the Taylor expansion. It buys quadratic convergence.
Monotone lifting on $[0, 1]$ with dyadic coefficients Selects $c = \tfrac{1}{2}$ Forces $p(\sigma) \ge \sigma$ on $[0, 1]$ and lands $(a, b, c)$ on dyadic rationals.

The unique solution satisfying all three is $(a, b, c) = (2, -\tfrac{3}{2}, \tfrac{1}{2})$, giving the polynomial

$$p(\sigma) = 2\sigma - \tfrac{3}{2}\sigma^3 + \tfrac{1}{2}\sigma^5.$$
(★)

This isn't a tuned hyperparameter. It's the answer to a well-posed question. The contrast with Stage 1 is what makes the whole design click. Stage 1's $(3.4445, -4.7750, 2.0315)$ comes out of a minimax fit, irrational, optimized for squeeze ratio across a wide input interval. Stage 2's $(2, -\tfrac{3}{2}, \tfrac{1}{2})$ comes out of algebra: three conditions, dyadic rationals, exact. Together they're a clean instance of an old engineering pattern: approximate coarsely, then correct exactly.

Three conditions, one triple. $(2, -\tfrac{3}{2}, \tfrac{1}{2})$ wasn't picked. It was forced.