Formula & theorem reference

Take an acute angle \theta, so its terminal side crosses the unit circle at a point (x,y) in the first quadrant. Drop a perpendicular segment straight down from (x,y) to the x-axis. That segment, the x-axis, and the terminal side form a right triangle with:

• hypotenuse = 1 (the radius of the unit circle),
• the side adjacent to \theta running along the x-axis, with length x,
• the side opposite \theta standing vertically, with length y.

By the right-triangle definition already in hand, \cos\theta = \frac{\text{adjacent}}{\text{hypotenuse}} = \frac{x}{1} = x, \qquad \sin\theta = \frac{\text{opposite}}{\text{hypotenuse}} = \frac{y}{1} = y. So for every acute angle, the coordinate definition and the right-triangle definition hand back the exact same numbers — the unit circle isn't replacing the old definition, it's the same idea with the hypotenuse fixed at 1 so "adjacent" and "opposite" collapse into plain x and y.

For \theta in any other quadrant, drop the same perpendicular from (x,y) to the x-axis and let \theta' be the resulting reference angle — the acute angle between the terminal side and the x-axis. That perpendicular still forms a right triangle with hypotenuse 1, legs |x| and |y|, and angle \theta', so the acute case already proven gives \cos\theta' = |x| and \sin\theta' = |y|. The sign is supplied automatically, not separately argued: x and y are the point's actual (signed) coordinates, negative exactly in the quadrants where the terminal side puts them, so \cos\theta = x and \sin\theta = y hold with their correct sign built in — nothing about the coordinate definition ever stops applying, only the right-triangle language ("opposite," "adjacent") stops making sense once \theta passes 90°.

Place triangle ABC with angle C at vertex C, flanked by sides CA = b and CB = a; the side opposite C is AB = c. Drop an altitude from vertex B down to line CA, meeting it at point D, with height h = BD.

Triangle BDC is a right triangle with hypotenuse CB = a and angle C at vertex C, so CD = a\cos C, \qquad h = a\sin C. (If C is obtuse, \cos C is negative and D falls past C, outside segment CA — the signed length CD = a\cos C still comes out correct, so nothing below needs to change for that case.)

That leaves AD = CA - CD = b - a\cos C. Triangle ABD has its right angle at D and hypotenuse AB = c, so the Pythagorean theorem gives c^2 = AD^2 + h^2 = (b - a\cos C)^2 + (a\sin C)^2. Expand and group the a^2 terms together: c^2 = b^2 - 2ab\cos C + a^2\cos^2 C + a^2\sin^2 C = b^2 - 2ab\cos C + a^2\left(\cos^2 C + \sin^2 C\right). The Pythagorean identity collapses \cos^2 C + \sin^2 C to 1: c^2 = a^2 + b^2 - 2ab\cos C. Relabeling the same argument around the other two vertices gives the A- and B-rotations.

Sanity check: set C = 90°. Then \cos C = 0, the whole correction term vanishes, and what's left is c^2 = a^2 + b^2 — the Pythagorean theorem, exactly as it must be, since a right triangle is just the C = 90° special case of every triangle.

Write the vectors in components, \mathbf{u} = (u_1,\dots,u_n), \mathbf{v} = (v_1,\dots,v_n), \mathbf{w} = (w_1,\dots,w_n), so that the dot product is the entrywise sum \mathbf{u}\cdot\mathbf{v} = \sum_{i=1}^n u_i v_i. Each property is now a statement about ordinary real numbers, one index at a time.

Commutativity. Multiplication of real numbers commutes, so u_i v_i = v_i u_i for every i. Summing both sides over i gives \mathbf{u}\cdot\mathbf{v} = \sum_i u_i v_i = \sum_i v_i u_i = \mathbf{v}\cdot\mathbf{u}.

Distributivity. The i-th entry of \mathbf{v}+\mathbf{w} is v_i + w_i, so \mathbf{u}\cdot(\mathbf{v}+\mathbf{w}) = \sum_{i=1}^n u_i(v_i+w_i) = \sum_{i=1}^n \left(u_i v_i + u_i w_i\right) = \sum_{i=1}^n u_i v_i + \sum_{i=1}^n u_i w_i = \mathbf{u}\cdot\mathbf{v} + \mathbf{u}\cdot\mathbf{w}, where the middle step splits the sum using ordinary distributivity and additivity of \sum.

Scalar compatibility. The i-th entry of c\mathbf{u} is cu_i, so (c\mathbf{u})\cdot\mathbf{v} = \sum_{i=1}^n (cu_i)v_i = c\sum_{i=1}^n u_i v_i = c(\mathbf{u}\cdot\mathbf{v}), pulling the constant c out of every term and then out of the sum. The same computation with \mathbf{v} scaled instead of \mathbf{u} gives \mathbf{u}\cdot(c\mathbf{v}) = c(\mathbf{u}\cdot\mathbf{v}), and combining the two chains of equalities yields all three expressions equal.

Positive definiteness. Taking \mathbf{v} = \mathbf{u} in the definition gives \mathbf{u}\cdot\mathbf{u} = \sum_{i=1}^n u_i^2, a sum of squares of real numbers, so \mathbf{u}\cdot\mathbf{u} \ge 0. If \mathbf{u} = \mathbf{0} then every u_i = 0 and the sum is 0. Conversely, if \sum_{i=1}^n u_i^2 = 0, then since each term u_i^2 \ge 0, every term must itself be 0 (a sum of nonnegative numbers vanishes only when each summand does), forcing u_i = 0 for every i and hence \mathbf{u} = \mathbf{0}.

This last property is what lets \lVert\mathbf{u}\rVert = \sqrt{\mathbf{u}\cdot\mathbf{u}} serve as a genuine length: the quantity under the root is never negative, and it is zero only for the zero vector itself.

If \mathbf{v} = \mathbf{0} both sides are 0 and the inequality holds trivially, so assume \mathbf{v} \neq \mathbf{0}.

Build a function of a single real variable t by looking at the squared length of \mathbf{u} - t\mathbf{v}. A squared length can never be negative, so for every t, 0 \le \lVert\mathbf{u} - t\mathbf{v}\rVert^2 = (\mathbf{u} - t\mathbf{v})\cdot(\mathbf{u} - t\mathbf{v}).

Expand the dot product using bilinearity: 0 \le \mathbf{u}\cdot\mathbf{u} - 2t\,(\mathbf{u}\cdot\mathbf{v}) + t^2\,(\mathbf{v}\cdot\mathbf{v}) = \lVert\mathbf{v}\rVert^2\,t^2 - 2(\mathbf{u}\cdot\mathbf{v})\,t + \lVert\mathbf{u}\rVert^2.

This is a quadratic in t with positive leading coefficient \lVert\mathbf{v}\rVert^2 (since \mathbf{v} \neq \mathbf{0}), and it never dips below 0 — so it can have at most one real root, meaning its discriminant cannot be positive: \left(-2(\mathbf{u}\cdot\mathbf{v})\right)^2 - 4\lVert\mathbf{v}\rVert^2\lVert\mathbf{u}\rVert^2 \le 0.

Simplify: 4(\mathbf{u}\cdot\mathbf{v})^2 \le 4\lVert\mathbf{u}\rVert^2\lVert\mathbf{v}\rVert^2 \quad\Longrightarrow\quad (\mathbf{u}\cdot\mathbf{v})^2 \le \lVert\mathbf{u}\rVert^2\lVert\mathbf{v}\rVert^2.

Taking square roots of both (nonnegative) sides gives |\mathbf{u}\cdot\mathbf{v}| \le \lVert\mathbf{u}\rVert\,\lVert\mathbf{v}\rVert.

Equality case. The discriminant is exactly 0 precisely when the quadratic touches 0 at its unique root t_0, i.e. when \lVert\mathbf{u} - t_0\mathbf{v}\rVert^2 = 0 for some t_0 — which forces \mathbf{u} - t_0\mathbf{v} = \mathbf{0}, so \mathbf{u} = t_0\mathbf{v}. Equality in the inequality therefore holds exactly when \mathbf{u} is a scalar multiple of \mathbf{v} (covering \mathbf{u} = \mathbf{0} as the t_0 = 0 case).

Sanity check: take \mathbf{u} = \mathbf{v}. Then the left side is |\mathbf{u}\cdot\mathbf{u}| = \lVert\mathbf{u}\rVert^2 and the right side is \lVert\mathbf{u}\rVert\,\lVert\mathbf{u}\rVert = \lVert\mathbf{u}\rVert^2 — equality, exactly as the parallel case demands since a vector is trivially a scalar multiple of itself.

Both sides are nonnegative, so it's enough to compare their squares. Expand the left side using the norm-squared identity \lVert \mathbf{w} \rVert^2 = \mathbf{w} \cdot \mathbf{w} and bilinearity of the dot product: \lVert \mathbf{u} + \mathbf{v} \rVert^2 = (\mathbf{u} + \mathbf{v})\cdot(\mathbf{u} + \mathbf{v}) = \mathbf{u}\cdot\mathbf{u} + 2(\mathbf{u}\cdot\mathbf{v}) + \mathbf{v}\cdot\mathbf{v} = \lVert \mathbf{u} \rVert^2 + 2(\mathbf{u}\cdot\mathbf{v}) + \lVert \mathbf{v} \rVert^2. Meanwhile the square of the right side is \left(\lVert \mathbf{u} \rVert + \lVert \mathbf{v} \rVert\right)^2 = \lVert \mathbf{u} \rVert^2 + 2\lVert \mathbf{u} \rVert\lVert \mathbf{v} \rVert + \lVert \mathbf{v} \rVert^2. The two expressions share the terms \lVert \mathbf{u} \rVert^2 and \lVert \mathbf{v} \rVert^2, so the whole inequality collapses to a claim about the cross terms alone: \mathbf{u}\cdot\mathbf{v} \le \lVert \mathbf{u} \rVert\lVert \mathbf{v} \rVert. This is the Cauchy–Schwarz inequality, and it's what actually needs proving.

Proving the cross-term bound. If \mathbf{v} = \mathbf{0} the claim reads 0 \le 0 and there's nothing to check, so assume \mathbf{v} \ne \mathbf{0}. For any real number t, the vector \mathbf{u} - t\mathbf{v} has nonnegative squared length, so the scalar function f(t) = \lVert \mathbf{u} - t\mathbf{v} \rVert^2 = (\mathbf{u} - t\mathbf{v})\cdot(\mathbf{u} - t\mathbf{v}) = \lVert \mathbf{v} \rVert^2 t^2 - 2(\mathbf{u}\cdot\mathbf{v})\,t + \lVert \mathbf{u} \rVert^2 satisfies f(t) \ge 0 for every t. This is a quadratic in t with positive leading coefficient \lVert \mathbf{v} \rVert^2 (since \mathbf{v} \ne \mathbf{0}), and a quadratic with positive leading coefficient that never dips below 0 can have at most one real root — its discriminant cannot be positive. Writing the discriminant out, \left(-2(\mathbf{u}\cdot\mathbf{v})\right)^2 - 4\lVert \mathbf{v} \rVert^2\lVert \mathbf{u} \rVert^2 \le 0, which simplifies to 4(\mathbf{u}\cdot\mathbf{v})^2 \le 4\lVert \mathbf{u} \rVert^2\lVert \mathbf{v} \rVert^2, hence (\mathbf{u}\cdot\mathbf{v})^2 \le \lVert \mathbf{u} \rVert^2\lVert \mathbf{v} \rVert^2. Taking square roots gives |\mathbf{u}\cdot\mathbf{v}| \le \lVert \mathbf{u} \rVert\lVert \mathbf{v} \rVert, and in particular \mathbf{u}\cdot\mathbf{v} \le \lVert \mathbf{u} \rVert\lVert \mathbf{v} \rVert, exactly the bound the first part of the argument was waiting on.

Closing the loop. Substituting this bound back into the expanded square of \lVert \mathbf{u} + \mathbf{v} \rVert gives \lVert \mathbf{u} + \mathbf{v} \rVert^2 = \lVert \mathbf{u} \rVert^2 + 2(\mathbf{u}\cdot\mathbf{v}) + \lVert \mathbf{v} \rVert^2 \le \lVert \mathbf{u} \rVert^2 + 2\lVert \mathbf{u} \rVert\lVert \mathbf{v} \rVert + \lVert \mathbf{v} \rVert^2 = \left(\lVert \mathbf{u} \rVert + \lVert \mathbf{v} \rVert\right)^2. Both \lVert \mathbf{u} + \mathbf{v} \rVert and \lVert \mathbf{u} \rVert + \lVert \mathbf{v} \rVert are nonnegative, so taking square roots preserves the inequality: \lVert \mathbf{u} + \mathbf{v} \rVert \le \lVert \mathbf{u} \rVert + \lVert \mathbf{v} \rVert.

Sanity check: if \mathbf{u} and \mathbf{v} point in exactly the same direction, say \mathbf{v} = c\mathbf{u} for some c > 0, then \mathbf{u}\cdot\mathbf{v} = c\lVert \mathbf{u} \rVert^2 = \lVert \mathbf{u} \rVert\lVert \mathbf{v} \rVert exactly, so every \le above becomes = and the two sides of the triangle inequality coincide — the shortcut and the detour are the same walk. Any other relative direction only shortens the straight path further, which is exactly the inequality's content: going straight from tail to tip never takes longer than going by way of an intermediate stop.

The whole argument comes from one identity: for any vectors \mathbf{u}, \mathbf{v}, \lVert\mathbf{u}+\mathbf{v}\rVert^2 = (\mathbf{u}+\mathbf{v})\cdot(\mathbf{u}+\mathbf{v}). Expand the right side using the distributive and commutative properties of the dot product: (\mathbf{u}+\mathbf{v})\cdot(\mathbf{u}+\mathbf{v}) = \mathbf{u}\cdot\mathbf{u} + \mathbf{u}\cdot\mathbf{v} + \mathbf{v}\cdot\mathbf{u} + \mathbf{v}\cdot\mathbf{v} = \mathbf{u}\cdot\mathbf{u} + 2(\mathbf{u}\cdot\mathbf{v}) + \mathbf{v}\cdot\mathbf{v}. Since \mathbf{u}\cdot\mathbf{u} = \lVert\mathbf{u}\rVert^2 and \mathbf{v}\cdot\mathbf{v} = \lVert\mathbf{v}\rVert^2, this becomes the vector law of cosines: \lVert\mathbf{u}+\mathbf{v}\rVert^2 = \lVert\mathbf{u}\rVert^2 + \lVert\mathbf{v}\rVert^2 + 2(\mathbf{u}\cdot\mathbf{v}).

That last equation holds unconditionally — no assumption about the angle between \mathbf{u} and \mathbf{v} was ever used. It says the gap between \lVert\mathbf{u}+\mathbf{v}\rVert^2 and \lVert\mathbf{u}\rVert^2 + \lVert\mathbf{v}\rVert^2 is exactly the term 2(\mathbf{u}\cdot\mathbf{v}), nothing more and nothing less. So the two sides agree precisely when that term vanishes: \lVert\mathbf{u}+\mathbf{v}\rVert^2 = \lVert\mathbf{u}\rVert^2 + \lVert\mathbf{v}\rVert^2 \quad\Longleftrightarrow\quad \mathbf{u}\cdot\mathbf{v} = 0 \quad\Longleftrightarrow\quad \mathbf{u} \perp \mathbf{v}, which is the two-way statement above — an "only if" as much as an "if," since every step from the expansion onward is a chain of equivalences, not just implications.

Sanity check: picture \mathbf{u} and \mathbf{v} as two legs of a right angle at the origin, so \mathbf{u}+\mathbf{v} is the vector running from the tip of one leg to the tip of the other along the hypotenuse of the triangle they span. The theorem then reads \lVert\mathbf{u}+\mathbf{v}\rVert^2 = \lVert\mathbf{u}\rVert^2 + \lVert\mathbf{v}\rVert^2 — hypotenuse squared equals the sum of the two legs squared, the ordinary Pythagorean theorem, now recovered as the \mathbf{u}\cdot\mathbf{v}=0 case of a statement that makes sense for vectors in any \mathbb{R}^n, not just the plane.

Write \mathbf{u} = (u_1, u_2, u_3) and \mathbf{v} = (v_1, v_2, v_3), and recall the defining formula \mathbf{u}\times\mathbf{v} = (u_2 v_3 - u_3 v_2,\ u_3 v_1 - u_1 v_3,\ u_1 v_2 - u_2 v_1).

Anticommutativity. Swapping \mathbf{u} and \mathbf{v} swaps every subscript pair inside each component: the first entry becomes v_2 u_3 - v_3 u_2, which is exactly -(u_2 v_3 - u_3 v_2). The same sign flip happens in the other two entries, so \mathbf{v}\times\mathbf{u} = -(\mathbf{u}\times\mathbf{v}).

Bilinearity. Each component of \mathbf{u}\times\mathbf{v} is a sum of products of coordinates, and coordinates of a sum add coordinatewise ((\mathbf{v}+\mathbf{w})_i = v_i + w_i), so plugging \mathbf{v}+\mathbf{w} into the formula above and distributing the products component by component recovers \mathbf{u}\times\mathbf{v} + \mathbf{u}\times\mathbf{w} term for term. Scaling one factor by c multiplies every product in every component by c, which factors straight back out, giving (c\mathbf{u})\times\mathbf{v} = c(\mathbf{u}\times\mathbf{v}).

Orthogonality. Compute the dot product directly: (\mathbf{u}\times\mathbf{v})\cdot\mathbf{u} = u_1(u_2 v_3 - u_3 v_2) + u_2(u_3 v_1 - u_1 v_3) + u_3(u_1 v_2 - u_2 v_1). Expanding gives six terms, u_1 u_2 v_3 - u_1 u_3 v_2 + u_2 u_3 v_1 - u_1 u_2 v_3 + u_1 u_3 v_2 - u_2 u_3 v_1, and each term appears twice with opposite sign — u_1u_2v_3 cancels -u_1u_2v_3, and likewise for the other two pairs — so the sum collapses to 0. The same cancellation pattern appears when dotting with \mathbf{v} instead, so (\mathbf{u}\times\mathbf{v})\cdot\mathbf{v} = 0 as well. A vector orthogonal to both \mathbf{u} and \mathbf{v} is normal to any plane containing them, which is the geometric content of the cross product.

Magnitude. Lagrange's identity states \lVert\mathbf{u}\times\mathbf{v}\rVert^2 = \lVert\mathbf{u}\rVert^2\lVert\mathbf{v}\rVert^2 - (\mathbf{u}\cdot\mathbf{v})^2, which a direct (if tedious) expansion of both sides into coordinates confirms term for term. Substituting \mathbf{u}\cdot\mathbf{v} = \lVert\mathbf{u}\rVert\lVert\mathbf{v}\rVert\cos\theta turns the right side into \lVert\mathbf{u}\rVert^2\lVert\mathbf{v}\rVert^2(1 - \cos^2\theta) = \lVert\mathbf{u}\rVert^2\lVert\mathbf{v}\rVert^2\sin^2\theta by the Pythagorean identity, and since 0 \le \theta \le \pi makes \sin\theta \ge 0, taking square roots gives \lVert\mathbf{u}\times\mathbf{v}\rVert = \lVert\mathbf{u}\rVert\lVert\mathbf{v}\rVert\sin\theta.

Parallelism test. If \mathbf{u} = \mathbf{0}, \mathbf{v} = \mathbf{0}, or \mathbf{v} = k\mathbf{u} for some scalar k, then \theta = 0 or \theta = \pi (or the magnitude formula's factor \lVert\mathbf{u}\rVert or \lVert\mathbf{v}\rVert is already 0), so \sin\theta = 0 and the magnitude formula forces \lVert\mathbf{u}\times\mathbf{v}\rVert = 0, hence \mathbf{u}\times\mathbf{v} = \mathbf{0}. Conversely, if neither vector is \mathbf{0} and they are not parallel, then 0 < \theta < \pi gives \sin\theta > 0, so the magnitude formula forces \lVert\mathbf{u}\times\mathbf{v}\rVert > 0 and \mathbf{u}\times\mathbf{v} cannot be \mathbf{0}.

Let P_0=(x_0,y_0,z_0) be any point that IS on the plane, so ax_0+by_0+cz_0=d. The distance D is the length of whatever piece of \overrightarrow{P_0P_1} points straight off the plane — exactly the absolute value of the scalar projection of \overrightarrow{P_0P_1} onto the unit normal \hat{\mathbf{n}}=\mathbf{n}/|\mathbf{n}|:

Substituting ax_0+by_0+cz_0=d gives D=\dfrac{|ax_1+by_1+cz_1-d|}{|\mathbf{n}|}, which no longer references P_0 at all — any point on the plane gives the same answer, as it must.

Build the parallelepiped one edge direction at a time and track how each step changes both sides of the claim.

One edge. A single vector \mathbf{u} spans a degenerate "box" — a line segment — with no volume, and a 1\times1 determinant argument doesn't apply here since we're in \mathbb{R}^3; skip straight to two edges, where the geometry first becomes 2-dimensional.

Two edges. \mathbf{u} and \mathbf{v} span a parallelogram in some plane through the origin. Decompose \mathbf{v} = \mathbf{v}_\parallel + \mathbf{v}_\perp, where \mathbf{v}_\parallel is the component of \mathbf{v} along \mathbf{u} and \mathbf{v}_\perp is the component orthogonal to \mathbf{u} within that plane. The parallelogram's area is base times height, A = \lVert\mathbf{u}\rVert \lVert\mathbf{v}_\perp\rVert, since \mathbf{v}_\parallel contributes no height above the line through \mathbf{u}.

Third edge. Now bring in \mathbf{w}, which is not required to lie in that plane. Split \mathbf{w} = \mathbf{w}_\parallel + \mathbf{w}_\perp relative to the plane spanned by \mathbf{u}, \mathbf{v}: \mathbf{w}_\perp is the component of \mathbf{w} orthogonal to that plane, and \mathbf{w}_\parallel lies inside it. Sliding \mathbf{w} by any multiple of \mathbf{u} or \mathbf{v} (i.e. replacing \mathbf{w} by \mathbf{w} - \mathbf{w}_\parallel = \mathbf{w}_\perp) shears the parallelepiped without changing its volume, because it slides the top face sideways parallel to the base while leaving the base and the perpendicular height unchanged. So V = (\text{area of base}) \times (\text{height}) = \left(\lVert\mathbf{u}\rVert \lVert\mathbf{v}_\perp\rVert\right) \lVert\mathbf{w}_\perp\rVert.

Matching the determinant to the same shear. Column operations of the form "add a multiple of one column to another" leave a determinant unchanged, exactly the algebraic mirror of the geometric shear above. Replacing the third column \mathbf{w} by \mathbf{w}_\perp = \mathbf{w} - \mathbf{w}_\parallel is such an operation (it subtracts the appropriate combination of the \mathbf{u}- and \mathbf{v}-columns), so \det\begin{bmatrix} \mathbf{u} & \mathbf{v} & \mathbf{w} \end{bmatrix} = \det\begin{bmatrix} \mathbf{u} & \mathbf{v} & \mathbf{w}_\perp \end{bmatrix}. Likewise, replacing \mathbf{v} by \mathbf{v}_\perp inside that determinant is a column operation that subtracts a multiple of the \mathbf{u}-column, so \det\begin{bmatrix} \mathbf{u} & \mathbf{v} & \mathbf{w}_\perp \end{bmatrix} = \det\begin{bmatrix} \mathbf{u} & \mathbf{v}_\perp & \mathbf{w}_\perp \end{bmatrix}. The three columns \mathbf{u}, \mathbf{v}_\perp, \mathbf{w}_\perp are now mutually orthogonal: \mathbf{v}_\perp \perp \mathbf{u} by construction, and \mathbf{w}_\perp is orthogonal to the entire plane spanned by \mathbf{u} and \mathbf{v} (hence to both \mathbf{u} and \mathbf{v}_\perp). Rescaling each of the three columns to a unit vector turns the matrix into an orthogonal matrix Q, and \det\begin{bmatrix} \mathbf{u} & \mathbf{v}_\perp & \mathbf{w}_\perp \end{bmatrix} = \lVert\mathbf{u}\rVert \lVert\mathbf{v}_\perp\rVert \lVert\mathbf{w}_\perp\rVert \det(Q), and \det(Q) = \pm 1 for any orthogonal matrix. Taking absolute values, \left|\det\begin{bmatrix} \mathbf{u} & \mathbf{v} & \mathbf{w} \end{bmatrix}\right| = \lVert\mathbf{u}\rVert \lVert\mathbf{v}_\perp\rVert \lVert\mathbf{w}_\perp\rVert, which is exactly the volume V computed geometrically above.

Sanity check: take \mathbf{u} = (1,0,0), \mathbf{v} = (0,1,0), \mathbf{w} = (0,0,1) — the unit cube. Every cross-term vanishes, the matrix is already the identity, and |\det I| = 1, matching the unit cube's volume of 1 exactly.

Left-multiplication reproduces the operation. Write I_n's rows as \mathbf{e}_1^{\mathsf T},\dots,\mathbf{e}_n^{\mathsf T} (each \mathbf{e}_i the i-th standard basis vector), so row i of any product CA is \mathbf{e}_i^{\mathsf T}A, a linear combination of the rows of A with coefficients read off from row i of C. This holds for every row of every product CA, so it holds in particular when C = E. Check each operation in turn:

Swap. E is I_n with rows p and q exchanged, so row p of E is \mathbf{e}_q^{\mathsf T} and row q of E is \mathbf{e}_p^{\mathsf T}, with every other row i still \mathbf{e}_i^{\mathsf T}. Row p of EA is then \mathbf{e}_q^{\mathsf T}A — row q of A — and row q of EA is row p of A, while every untouched row of EA equals the matching row of A. That's exactly rows p and q of A swapped.

Scale. E is I_n with row p replaced by k\mathbf{e}_p^{\mathsf T} for some k \ne 0. Row p of EA is k\mathbf{e}_p^{\mathsf T}A = k times row p of A, and every other row of E is still \mathbf{e}_i^{\mathsf T}, so every other row of EA matches A unchanged.

Row combination. E is I_n with row p replaced by \mathbf{e}_p^{\mathsf T} + k\mathbf{e}_q^{\mathsf T} for some k and some q \ne p. Row p of EA is then \mathbf{e}_p^{\mathsf T}A + k\mathbf{e}_q^{\mathsf T}A, which is row p of A plus k times row q of A — every other row of E, and hence of EA, is untouched.

In all three cases EA is precisely A with the one named row operation applied, which is the claim.

Invertibility. Each elementary row operation has an inverse operation of the same three kinds: undo a swap of rows p,q by swapping p,q again; undo scaling row p by k \ne 0 by scaling row p by 1/k; undo adding k times row q to row p by adding -k times row q to row p. Let F be the elementary matrix for that inverse operation. Applying E's operation to I_n produces E by definition, and applying F's operation to E undoes it, landing back at I_n; by the left-multiplication fact just proved (with A = E), that composite is FE. So FE = I_n, and running the two operations in the opposite order gives EF = I_n by the identical argument. Hence E is invertible with E^{-1} = F, and F is elementary by construction.

Sanity check: take n=2 and E = \begin{pmatrix}1&0\\-3&1\end{pmatrix} — the operation "add -3 times row 1 to row 2." Its inverse operation is "add +3 times row 1 to row 2," giving F=\begin{pmatrix}1&0\\3&1\end{pmatrix}, and indeed FE=\begin{pmatrix}1&0\\3&1\end{pmatrix}\begin{pmatrix}1&0\\-3&1\end{pmatrix} =\begin{pmatrix}1&0\\3-3&1\end{pmatrix}=I_2, confirming E^{-1}=F.

Track the elimination as a sequence of matrix multiplications rather than as row-by-row bookkeeping. Each row replacement step — "add m times row j to row i" — is itself left-multiplication by an elementary matrix E: the identity matrix with m dropped into the (i,j) slot. If reducing A to U takes k such steps with elementary matrices E_1, E_2, \dots, E_k applied in that order, then

Every E_t here is unit lower triangular (1's on the diagonal, a single nonzero entry below it, from a row replacement that only ever adds a multiple of a higher row to a lower one — no swaps means row i only ever absorbs rows above it). The set of unit lower triangular matrices is closed under multiplication and under inversion, so the product M = E_k \cdots E_2 E_1 is itself unit lower triangular, and so is its inverse L = M^{-1}. Solving MA = U for A gives

It remains to identify the entries of L. Each E_t^{-1} undoes its row operation by flipping the sign of its one off-diagonal entry: if E_t adds m times row j to row i, then E_t^{-1} subtracts it back out, so E_t^{-1} carries -m in slot (i,j). Write \ell_{ij} = -m for that entry, i.e. \ell_{ij} is the negative of the multiplier used to clear (i,j) — equivalently, the multiplier used to reconstruct it from U's reduction of A. Because no swaps occur, the elimination order never revisits a cleared entry, so the individual E_t^{-1} matrices combine without their off-diagonal entries interfering with one another, and

comes out unit lower triangular with exactly \ell_{ij} sitting in each below-diagonal slot (i,j) — precisely the multiplier record from the forward elimination. This is why L can be written down directly from the row reduction as it happens, with no matrix inversion required in practice.

Sanity check: if A is already upper triangular, elimination needs no steps at all, so L = I and U = A — consistent with A = IA.

Start from the Leibniz definition of the determinant as a sum over every permutation \sigma of \{1, \dots, n\}: \det(A) = \sum_{\sigma \in S_n} \operatorname{sgn}(\sigma) \prod_{k=1}^n a_{k,\sigma(k)}. Fix a row i and sort this sum by which column each permutation sends row i to. Every \sigma \in S_n satisfies \sigma(i) = j for exactly one j \in \{1, \dots, n\}, so the sum splits cleanly into n disjoint groups: \det(A) = \sum_{j=1}^n \sum_{\substack{\sigma \in S_n \\ \sigma(i)=j}} \operatorname{sgn}(\sigma) \prod_{k=1}^n a_{k,\sigma(k)}. In the inner sum, every term already carries the factor a_{i,\sigma(i)} = a_{ij}, which does not depend on \sigma beyond fixing \sigma(i) = j, so it factors out: \det(A) = \sum_{j=1}^n a_{ij} \left( \sum_{\substack{\sigma \in S_n \\ \sigma(i)=j}} \operatorname{sgn}(\sigma) \prod_{k \neq i} a_{k,\sigma(k)} \right). The claim is that the bracketed sum is exactly the cofactor C_{ij}. Each permutation \sigma with \sigma(i) = j restricts to a bijection between the remaining n-1 row indices \{1,\dots,n\}\setminus\{i\} and the remaining n-1 column indices \{1,\dots,n\}\setminus\{j\} — call this restricted map \tau. Relabeling those two index sets in increasing order identifies \tau with a permutation \bar\tau of \{1, \dots, n-1\}, and the product \prod_{k \neq i} a_{k,\sigma(k)} becomes exactly the corresponding product of entries of the submatrix obtained by deleting row i and column j — the submatrix whose determinant is M_{ij}. Summing \operatorname{sgn}(\bar\tau) \prod \bar\tau\text{-entries} over all such \bar\tau gives M_{ij} itself, provided the signs match up: that \operatorname{sgn}(\sigma) = (-1)^{i+j} \operatorname{sgn}(\bar\tau).

To see the sign correction, compare \sigma to the permutation \bar\tau extended back to \{1,\dots,n\} by sending i \mapsto j in the "obvious" position. Moving row i up to row 1 by swapping it past the i-1 rows above it takes i - 1 transpositions; moving column j left to column 1 by the same process takes j - 1 transpositions. Each such adjacent swap flips the sign of a permutation by -1, so relating \sigma's sign to \bar\tau's sign (computed on the un-shifted remaining indices) costs a factor of (-1)^{(i-1)+(j-1)} = (-1)^{i+j}. Hence \sum_{\substack{\sigma \in S_n \\ \sigma(i)=j}} \operatorname{sgn}(\sigma) \prod_{k \neq i} a_{k,\sigma(k)} = (-1)^{i+j} M_{ij} = C_{ij}. Substituting back gives the row-i expansion \det(A) = \sum_{j=1}^n a_{ij} C_{ij}. The column-j expansion follows by the identical argument applied to A^T, using \det(A) = \det(A^T) and the fact that transposing swaps the roles of rows and columns (and swaps M_{ij} \leftrightarrow M_{ji}, while (-1)^{i+j} = (-1)^{j+i} is unaffected).

Sanity check: for n = 2, expanding along row 1 gives \det(A) = a_{11}C_{11} + a_{12}C_{12} = a_{11}(a_{22}) - a_{12}(a_{21}) = a_{11}a_{22} - a_{12}a_{21} — the familiar 2\times 2 determinant formula, since each 1\times 1 minor is just the surviving entry itself.

Base case, n = 2. Write A = \begin{pmatrix}a&b\\c&d\end{pmatrix}, so \det A = ad - bc.

• Swap: \det\begin{pmatrix}c&d\\a&b\end{pmatrix} = cb - da = -(ad-bc).
• Scale row 1 by k: \det\begin{pmatrix}ka&kb\\c&d\end{pmatrix} = (ka)d - (kb)c = k(ad-bc).
• Replace row 2 with row 2 plus k times row 1: \det\begin{pmatrix}a&b\\c+ka&d+kb\end{pmatrix} = a(d+kb) - b(c+ka) = ad + kab - bc - kab = ad - bc.

Each identity matches the claim exactly, and scaling or replacing row 1 instead of row 2 (or vice versa) runs through the identical algebra with the roles of the rows exchanged. So all three properties hold for n=2.

Induction step. Suppose the three properties hold for every (n-1)\times(n-1) determinant, and let A be n\times n with n \ge 3. Every elementary row operation touches at most two rows of A, so with n \ge 3 rows on hand there is always some row index p that the operation leaves untouched. Expand \det A by cofactors along that row: \det A = \sum_{j=1}^{n} (-1)^{p+j}\, a_{pj}\, M_{pj}, where M_{pj} is the determinant of the (n-1)\times(n-1) minor formed by deleting row p and column j from A.

Because row p is untouched by the operation, the entries a_{pj} and the alternating sign pattern (-1)^{p+j} in this sum are identical for A and for B — only the minors can change. Deleting column j never removes either of the rows the operation acts on (that operation never involves row p, and row p is the one being deleted), so each minor M_{pj} simply carries the same elementary row operation, restricted to its own n-1 rows and n-1 columns. By the induction hypothesis:

• Swap: every minor M_{pj} negates, so every term (-1)^{p+j}a_{pj}M_{pj} negates, and summing negates \det A as a whole — \det B = -\det A.
• Scale: every minor M_{pj} scales by k (the scaled row survives deletion of column j intact, since only row p is deleted), so every term scales by k, and \det B = k\det A.
• Replacement: every minor M_{pj} is unchanged, so every term is unchanged, and \det B = \det A.

This closes the induction: the three properties hold for every n \ge 2.

Consequence. These are exactly the operations row reduction uses, so they justify computing a determinant by row-reducing A to an upper triangular matrix U and multiplying diagonal entries: each swap along the way contributes a factor -1, each row replacement contributes a factor 1, and \det A is recovered as \pm\det U (dividing out any scaling factors used). They also explain why a matrix with two equal or proportional rows has determinant 0: replacing one of those rows with itself minus the appropriate multiple of the other produces a zero row without changing the determinant by property 3, and a zero row forces determinant 0 (the k=0 case of property 2, since scaling any row by 0 both zeroes that row and scales \det A by 0).

Fix A and treat \det(AB) as a function of B alone, writing B = \begin{pmatrix} \mathbf{b}_1 & \mathbf{b}_2 & \cdots & \mathbf{b}_n \end{pmatrix} in terms of its columns, so that AB = \begin{pmatrix} A\mathbf{b}_1 & A\mathbf{b}_2 & \cdots & A\mathbf{b}_n \end{pmatrix}.

Case 1: B is singular. Then \det B = 0, and the columns of B are linearly dependent — some \mathbf{b}_k is a linear combination of the others. Because A acts linearly, A\mathbf{b}_k is that same combination of the columns A\mathbf{b}_1, \dots, A\mathbf{b}_n, so the columns of AB are dependent too. A matrix with dependent columns has determinant 0, so \det(AB) = 0 = \det(A) \cdot 0 = \det(A)\det(B). The identity holds in this case regardless of what A is.

Case 2: B is invertible. Then B is a product of elementary matrices, B = E_1 E_2 \cdots E_k, since row-reducing B to the identity records a sequence of elementary row operations, and undoing them (in reverse) rebuilds B from I one elementary factor at a time. Each elementary matrix E scales the determinant of whatever it multiplies by a fixed factor c_E depending only on which operation E performs — a row swap gives c_E = -1, scaling a row by r gives c_E = r, and adding a multiple of one row to another gives c_E = 1 — and in every case that same factor satisfies \det E = c_E \cdot \det I = c_E. So for any matrix M, \det(EM) = c_E \det M = (\det E)(\det M). Applying this one elementary factor at a time, \det(AB) = \det(A E_1 E_2 \cdots E_k) = (\det E_k)\det(AE_1\cdots E_{k-1}) = \cdots = \det(A)\det(E_1)\det(E_2)\cdots\det(E_k). Running the same one-factor-at-a-time collapse on B itself (with A replaced by I) gives \det(B) = \det(E_1)\det(E_2)\cdots\det(E_k), so the product of the \det(E_i)'s appearing above is exactly \det(B), giving \det(AB) = \det(A)\det(B).

Together the two cases cover every B, so the identity holds in general.

Sanity check: take A = B = I_n. Both sides read \det(I_n) = 1 \cdot 1, which checks out. Take B singular and A = I_n: the left side is \det(I_n B) = \det B = 0, matching \det(I_n)\det(B) = 1 \cdot 0.

Both directions lean on one algebraic fact borrowed from determinant theory: for any two n\times n matrices A and B, \det(AB) = (\det A)(\det B).

If A is invertible, then \det A \ne 0. By definition AA^{-1} = I_n. Taking the determinant of both sides and applying the product rule above, \det(A)\det\left(A^{-1}\right) = \det I_n = 1. A product of two real numbers can only equal 1 if neither factor is 0 (anything times 0 is 0, never 1), so \det A \ne 0 — and this same equation already isolates \det\left(A^{-1}\right) = \frac{1}{\det A}, which is the stated corollary.

If \det A \ne 0, then A is invertible. Row-reduce A to an echelon form U. Each elementary row operation used along the way changes \det A by a known nonzero factor — a row swap multiplies the determinant by -1, scaling a row by a nonzero constant c multiplies the determinant by c, and adding a multiple of one row to another leaves the determinant unchanged — so \det U is \det A times a nonzero constant. Since \det A \ne 0 by hypothesis, this forces \det U \ne 0 as well.

U is upper triangular, so \det U is exactly the product of its diagonal entries. For that product to be nonzero, every diagonal entry of U must be nonzero — none of them collapsed to a zero pivot during elimination. An n\times n triangular matrix with n nonzero diagonal entries has a pivot in every row and every column, so U row-reduces further, by scaling each row to make its pivot 1 and clearing every entry above each pivot, all the way down to I_n. Since A was row equivalent to U and U is row equivalent to I_n, A is row equivalent to I_n.

Every elementary row operation is itself invertible — undoing a swap is another swap, undoing a scaling by c\ne0 is scaling by 1/c, undoing a row addition is subtracting the same multiple back off — and each can be realized as left-multiplication by an elementary matrix E_k. Chaining the operations that carried A to I_n gives E_k \cdots E_1 A = I_n, so the product B = E_k\cdots E_1 satisfies BA = I_n. Since B is itself a product of invertible elementary matrices, B is invertible with inverse B^{-1} = E_1^{-1}\cdots E_k^{-1}; multiplying BA = I_n on the left by B^{-1} gives A = B^{-1}, so A is invertible with A^{-1} = B.

Sanity check: for the 1\times 1 case A = [a], \det A = a, and the theorem reduces to the familiar fact that a single number is invertible (has a reciprocal) exactly when it is nonzero — the determinant test is the n\times n generalization of that one observation.

Write A in terms of its columns, A = \begin{pmatrix} \mathbf{a}_1 & \mathbf{a}_2 & \cdots & \mathbf{a}_n \end{pmatrix}, and let \mathbf{x} = (x_1, \dots, x_n) be the solution of A\mathbf{x} = \mathbf{b}, which exists and is unique because A is invertible. Expanding the matrix product by columns, \mathbf{b} = A\mathbf{x} = x_1\mathbf{a}_1 + x_2\mathbf{a}_2 + \cdots + x_n\mathbf{a}_n.

Fix an index i and look at \det A_i(\mathbf{b}), the determinant with column i swapped out for \mathbf{b} and every other column left alone. Substitute the sum above for that one column: \det A_i(\mathbf{b}) = \det\begin{pmatrix} \mathbf{a}_1 & \cdots & \mathbf{a}_{i-1} & x_1\mathbf{a}_1 + \cdots + x_n\mathbf{a}_n & \mathbf{a}_{i+1} & \cdots & \mathbf{a}_n \end{pmatrix}.

The determinant is linear in each column separately, so this splits into a sum of n determinants, one for each term x_j\mathbf{a}_j: \det A_i(\mathbf{b}) = \sum_{j=1}^n x_j \det\begin{pmatrix} \mathbf{a}_1 & \cdots & \mathbf{a}_{i-1} & \mathbf{a}_j & \mathbf{a}_{i+1} & \cdots & \mathbf{a}_n \end{pmatrix}.

Look at each term on the right. When j \neq i, the matrix in that term has \mathbf{a}_j sitting in position i while \mathbf{a}_j also already occupies its own position j elsewhere among the untouched columns — the column list repeats \mathbf{a}_j twice. A matrix with two equal columns has determinant zero, so every j \neq i term vanishes. Only j = i survives, and in that term the matrix is just A itself (column i replaced by \mathbf{a}_i is no replacement at all), giving \det A_i(\mathbf{b}) = x_i \det A.

Since A is invertible, \det A \neq 0, so this can be solved for x_i: x_i = \frac{\det A_i(\mathbf{b})}{\det A}.

The index i was arbitrary, so this holds for every entry of \mathbf{x}.

Sanity check: for n = 2 with A = \begin{pmatrix} a & b \\ c & d \end{pmatrix} and \mathbf{b} = \begin{pmatrix} e \\ f \end{pmatrix}, the rule gives x_1 = \dfrac{ed - bf}{ad-bc} and x_2 = \dfrac{af - ec}{ad-bc} — exactly the formulas that come out of solving the 2\times2 system directly by elimination.

Fix a column index k and expand \det A by cofactors along column k: \det A = \sum_{i=1}^n a_{ik}C_{ik}. That's the ordinary cofactor-expansion formula — nothing new yet. The key step is what happens when the column being multiplied against the cofactors is the wrong one.

Claim: for j \neq k, \sum_{i=1}^n a_{ij}C_{ik} = 0. To see why, build a new matrix B by taking A and replacing column k with a copy of column j — so B has two identical columns (column j and column k both equal to A's column j), and every column of B other than column k agrees with A. A matrix with two identical columns has determinant 0: swapping those two columns leaves B unchanged, yet swapping any two columns of a matrix negates its determinant, so \det B = -\det B, forcing \det B = 0.

Now expand \det B by cofactors along column k. Deleting column k from B is exactly the same deletion as in A — column k is the one piece of B that differs from A, and it's the one being deleted — so every minor, and hence every cofactor C_{ik}, is identical for B and for A. The only thing that changed is the entries being multiplied against those cofactors: B's column k is A's column j. So 0 = \det B = \sum_{i=1}^n b_{ik}C_{ik} = \sum_{i=1}^n a_{ij}C_{ik}, which proves the claim.

Assembling the product A \cdot \operatorname{adj}(A). By definition the (j,k) entry of \operatorname{adj}(A) is C_{kj}, so the (j,k) entry of the product A\operatorname{adj}(A) is \sum_{i=1}^n a_{ji}\big[\operatorname{adj}(A)\big]_{ik} = \sum_{i=1}^n a_{ji}C_{ki}. Relabeling the summation index doesn't change which sum this is — it's exactly \sum_i a_{ji}C_{ki}, a sum of row-j entries against row-k cofactors. Cofactor expansion along a row obeys the identical two identities just proved for columns (the argument is the mirror image, using row swaps instead of column swaps): this sum equals \det A when j = k, and equals 0 when j \neq k. So A\operatorname{adj}(A) = (\det A)I_n, every off-diagonal entry vanishing by the claim and every diagonal entry equal to \det A by ordinary cofactor expansion along that row.

Since \det A \neq 0, divide both sides by \det A: A\left(\frac{1}{\det A}\operatorname{adj}(A)\right) = I_n. An n\times n matrix with a right inverse of matching size is invertible with that right inverse as its two-sided inverse, so A^{-1} = \frac{1}{\det A}\operatorname{adj}(A).

Sanity check (n=2): for A = \begin{pmatrix} a & b \\ c & d \end{pmatrix}, the cofactors are C_{11}=d, C_{12}=-c, C_{21}=-b, C_{22}=a, so \operatorname{adj}(A) = \begin{pmatrix} d & -b \\ -c & a \end{pmatrix} — the transpose of the cofactor matrix — and the formula reduces to A^{-1} = \frac{1}{ad-bc}\begin{pmatrix} d & -b \\ -c & a \end{pmatrix}, exactly the familiar 2\times 2 inverse formula.

(\implies) If W is a subspace, it is itself a vector space under the addition and scalar multiplication inherited from V, so both closure properties hold automatically — that direction is immediate from the definition of "subspace" and needs no further argument.

(\impliedby) This is the direction that does real work: assume W is nonempty and closed under addition and scalar multiplication, and check that W inherits every vector space axiom from V.

First, \mathbf{0} \in W. Since W is nonempty, pick any \mathbf{u} \in W. Closure under scalar multiplication with c = 0 gives 0 \cdot \mathbf{u} = \mathbf{0} \in W. So W contains the zero vector without that having to be assumed separately — it falls out of nonemptiness plus scalar closure.

Next, every \mathbf{u} \in W has its negative in W: closure under scalar multiplication with c = -1 gives (-1)\mathbf{u} = -\mathbf{u} \in W.

With \mathbf{0} \in W and negatives available inside W, the two closure hypotheses supply exactly the operations W needs: addition of any two elements of W stays in W, and scaling any element of W by any real number stays in W. So W is closed under both vector space operations, and it contains an additive identity and additive inverses drawn from that same closure.

Every remaining axiom — commutativity and associativity of addition, associativity of scalar multiplication, the distributive laws, and the identity 1 \cdot \mathbf{u} = \mathbf{u} — is an equation that holds for all vectors in V, by hypothesis that V is a vector space. Since every element of W is in particular an element of V, each of these equations holds for the elements of W too; W is not required to satisfy some different, weaker set of axioms, only the same ones V already satisfies for every vector it contains. So none of these need separate verification — they are inherited, not re-derived.

Having the two closure properties, an identity element, inverses, and every other axiom inherited for free, W satisfies the full definition of a vector space under the operations of V, restricted to W. That makes W a subspace of V.

Why the test stops at two conditions. A first guess might list three requirements — closure under addition, closure under scalar multiplication, and \mathbf{0} \in W — and that three-condition version is also correct and commonly used. The proof above shows why the third condition is redundant whenever W is already known to be nonempty: scaling any single element of W by 0 produces \mathbf{0}, so nonemptiness plus scalar closure already forces \mathbf{0} \in W. The nonemptiness hypothesis cannot be dropped, though — \varnothing satisfies both closure conditions vacuously (there are no pairs or elements to violate them) but is not a vector space, since it has no zero vector.

Sanity check: take V = \mathbb{R}^2 and W = \{(x, y) : y = x\}, the line through the origin. Adding two points on that line, (a, a) + (b, b) = (a+b, a+b), lands back on the line, and scaling, c(a,a) = (ca, ca), does too — both closure conditions hold, so the test confirms W is a subspace, matching the geometric fact that a line through the origin is a subspace of the plane. Shift the line off the origin to W' = \{(x, y) : y = x + 1\} and closure under scalar multiplication already fails: (0,1) \in W' but 2(0,1) = (0,2) \notin W', correctly flagging that a line missing the origin is not a subspace.

Write W = \operatorname{span}\{\mathbf{v}_1, \dots, \mathbf{v}_k\}, so W is exactly the set of all linear combinations c_1\mathbf{v}_1 + c_2\mathbf{v}_2 + \cdots + c_k\mathbf{v}_k with scalars c_1, \dots, c_k. Checking that W is a subspace of V means checking three things: \mathbf{0} \in W, W is closed under addition, and W is closed under scalar multiplication.

Contains the zero vector. Taking every coefficient to be 0 gives the linear combination 0\mathbf{v}_1 + 0\mathbf{v}_2 + \cdots + 0\mathbf{v}_k = \mathbf{0}, so \mathbf{0} \in W. In particular W is nonempty.

Closed under addition. Let \mathbf{u}, \mathbf{w} \in W. Then by definition of W there are scalars a_1, \dots, a_k and b_1, \dots, b_k with \mathbf{u} = a_1\mathbf{v}_1 + \cdots + a_k\mathbf{v}_k, \qquad \mathbf{w} = b_1\mathbf{v}_1 + \cdots + b_k\mathbf{v}_k. Adding these and grouping term by term, \mathbf{u} + \mathbf{w} = (a_1 + b_1)\mathbf{v}_1 + \cdots + (a_k + b_k)\mathbf{v}_k, which is again a linear combination of \mathbf{v}_1, \dots, \mathbf{v}_k — with coefficients a_i + b_i — so \mathbf{u} + \mathbf{w} \in W.

Closed under scalar multiplication. Let \mathbf{u} \in W as above and let t be any scalar. Then t\mathbf{u} = t(a_1\mathbf{v}_1 + \cdots + a_k\mathbf{v}_k) = (ta_1)\mathbf{v}_1 + \cdots + (ta_k)\mathbf{v}_k, again a linear combination of \mathbf{v}_1, \dots, \mathbf{v}_k, so t\mathbf{u} \in W.

All three conditions hold, so W is a subspace of V.

Smallest such subspace. Every \mathbf{v}_i lies in W, since \mathbf{v}_i is the linear combination with coefficient 1 on \mathbf{v}_i and 0 elsewhere. Now suppose U is any subspace of V with \mathbf{v}_1, \dots, \mathbf{v}_k \in U. Because U is closed under scalar multiplication, each c_i\mathbf{v}_i \in U, and because U is closed under addition, the sum c_1\mathbf{v}_1 + \cdots + c_k\mathbf{v}_k \in U for every choice of scalars. That sum is an arbitrary element of W, so W \subseteq U. Hence W sits inside every subspace containing the \mathbf{v}_i's, which is exactly what it means for W to be the smallest one.

Existence. Because \mathcal{B} is a basis for V, it spans V: \operatorname{span}(\mathcal{B}) = V. So \mathbf{x} \in V already guarantees some list of scalars c_1, \dots, c_n producing the combination above — that half costs nothing beyond the definition of span. The real content of the theorem is the word "exactly."

Uniqueness. Suppose \mathbf{x} has two such representations, \mathbf{x} = c_1\mathbf{b}_1 + \cdots + c_n\mathbf{b}_n \qquad\text{and}\qquad \mathbf{x} = d_1\mathbf{b}_1 + \cdots + d_n\mathbf{b}_n. Subtracting the second equation from the first leaves \mathbf{0} on the left: \mathbf{0} = (c_1 - d_1)\mathbf{b}_1 + (c_2 - d_2)\mathbf{b}_2 + \cdots + (c_n - d_n)\mathbf{b}_n. This is a linear combination of \mathbf{b}_1, \dots, \mathbf{b}_n equal to the zero vector. Since \mathcal{B} is a basis, its vectors are linearly independent, so the only linear combination of them that produces \mathbf{0} is the trivial one — every coefficient forced to 0: c_1 - d_1 = 0, \quad c_2 - d_2 = 0, \quad \dots, \quad c_n - d_n = 0. Hence c_i = d_i for every i = 1, \dots, n, so the two representations were the same list of scalars all along.

Why both properties are needed. Spanning alone (without independence) would let a redundant vector in \mathcal{B} soak up the same \mathbf{x} in more than one way — existence would survive but uniqueness would fail. Independence alone (without spanning) would keep every representation unique wherever one exists, but there could be vectors in V that \mathcal{B} never reaches at all. A basis is exactly the balance point: just enough vectors to reach all of V (spanning), and no more than that (independence) — which is precisely why "\mathbf{x} = c_1\mathbf{b}_1 + \cdots + c_n\mathbf{b}_n" pins down one and only one coefficient list, with no ambiguity about which scalars belong to it.

First claim. Suppose \mathbf{v}_k is a linear combination of the other vectors in S, say \mathbf{v}_k = c_1\mathbf{v}_1 + \cdots + c_{k-1}\mathbf{v}_{k-1} + c_{k+1}\mathbf{v}_{k+1} + \cdots + c_p\mathbf{v}_p. Write S' = S \setminus \{\mathbf{v}_k\}. Since S' \subseteq S, every vector reachable as a combination of S' is automatically reachable as a combination of S, so \operatorname{span}(S') \subseteq H.

For the reverse inclusion, take any \mathbf{y} \in H. Then \mathbf{y} is some combination \mathbf{y} = a_1\mathbf{v}_1 + \cdots + a_k\mathbf{v}_k + \cdots + a_p\mathbf{v}_p. Substitute the expression for \mathbf{v}_k above in place of \mathbf{v}_k and collect terms: each \mathbf{v}_i with i \neq k now carries the combined coefficient a_i + a_k c_i, and \mathbf{v}_k has vanished entirely. So \mathbf{y} is a linear combination of S' alone, giving H \subseteq \operatorname{span}(S'). The two inclusions together give \operatorname{span}(S') = H.

Second claim. If H \neq \{\mathbf{0}\}, then S contains at least one nonzero vector, so among all subsets of S that span H, choose one, call it B, with the fewest vectors possible (such a subset exists because S itself spans H and S is finite, so there are only finitely many subsets to compare). It remains to show B is linearly independent, which makes it a basis for H.

Suppose toward a contradiction that B is linearly dependent. Then some vector in B is a linear combination of the others in B — removing that vector, by the first claim applied to B in place of S, leaves a smaller subset of B (hence of S) that still spans H. This contradicts the minimality of B. So B is linearly independent, and since B already spans H, B is a basis for H.

Sanity check: with S = \{\mathbf{v}_1, \mathbf{v}_2\} where \mathbf{v}_2 = 3\mathbf{v}_1 and \mathbf{v}_1 \neq \mathbf{0}, the second vector is a multiple of the first, so dropping it leaves S' = \{\mathbf{v}_1\} with \operatorname{span}(S') = \operatorname{span}(S) exactly as the theorem predicts — and \{\mathbf{v}_1\} is itself the minimal, linearly independent spanning set, i.e. the basis the second claim guarantees.

Write S = \{v_1, \dots, v_n\}. Two facts about V are already on hand from the definition of dimension: every linearly independent subset of V can be grown, one vector at a time, into a basis of V; and every basis of V has exactly n vectors, since \dim V = n.

Independent \implies spanning. Suppose S is linearly independent but, toward a contradiction, does not span V. Then some w \in V lies outside \operatorname{span}(S). Adjoining w to S keeps the set independent: if c_1v_1 + \cdots + c_nv_n + c_{n+1}w = 0 with c_{n+1} \ne 0, solving for w would express it as a combination of v_1,\dots,v_n, placing w in \operatorname{span}(S) — which was excluded — so c_{n+1} = 0, and then the independence of S forces every remaining c_i = 0 as well. So S \cup \{w\} is an independent set of n+1 vectors. Extend it to a basis of V exactly as in the growing process above; that basis has at least n+1 vectors, contradicting \dim V = n. No such w exists, so S spans V, and an independent spanning set is by definition a basis.

Spanning \implies independent. Suppose S spans V but, toward a contradiction, is linearly dependent. Dependence means some v_k is a linear combination of the other n-1 vectors in S; deleting v_k leaves a set of n-1 vectors whose span is unchanged, since every occurrence of v_k in a combination can be rewritten using the rest. So V is spanned by only n - 1 vectors. A basis can always be extracted from a spanning set by discarding vectors that are combinations of the others already kept, so V would have a basis of at most n-1 vectors — contradicting \dim V = n. No such dependence exists, so S is independent, and again an independent spanning set is a basis.

Sanity check: in \mathbb{R}^2, where n = 2, the two vectors (1,1) and (1,-1) are independent (neither is a scalar multiple of the other), so the theorem guarantees — without any separate spanning check — that they already span all of \mathbb{R}^2.

Row reduce A to a reduced echelon form R. Row reduction is reversible (each step is undone by another row operation), so Ax = 0 and Rx = 0 have exactly the same solution set — hence \operatorname{Nul} A = \operatorname{Nul} R — and row reduction changes neither the number of pivot columns nor the column count, so \operatorname{rank}(A) = \operatorname{rank}(R) as well. It's therefore enough to prove the claim for R itself.

Say R has r pivot columns; by definition \operatorname{rank}(A) = r, and the remaining n - r columns are free. Solve Rx = 0 by the usual back-substitution: assign each free variable its own parameter t_1, \dots, t_{n-r}, and each pivot variable is then forced to some fixed linear combination of those parameters, read off from its row of R. Every solution vector x splits along these parameters as x = t_1 v_1 + t_2 v_2 + \cdots + t_{n-r} v_{n-r} for vectors v_1, \dots, v_{n-r} built entirely from the pivot rows' coefficients, so this set spans \operatorname{Nul} R.

The list is also independent: in the equation \sum_k t_k v_k = 0, look at entry j of v_k, where j is the position of the k-th free variable. By construction v_k has a 1 there and every other v_i (i \neq k) has a 0 there, since a free variable's own entry doesn't depend on the other free variables. So that entry of the sum reduces to t_k = 0, and this holds for every k. Hence v_1, \dots, v_{n-r} is a basis for \operatorname{Nul} R, giving \dim(\operatorname{Nul} R) = n - r.

Combining the two pieces, \operatorname{rank}(A) + \dim(\operatorname{Nul} A) = \operatorname{rank}(R) + \dim(\operatorname{Nul} R) = r + (n - r) = n.

Sanity check: an invertible n \times n matrix has \operatorname{Nul} A = \{0\}, so the theorem forces \operatorname{rank}(A) = n — every column pivots, exactly as invertibility demands.

Let k = \dim\left(\operatorname{Nul}A\right) and pick a basis \{\mathbf{v}_1,\dots,\mathbf{v}_k\} for \operatorname{Nul}A, the set of vectors in \mathbb{R}^n that A sends to \mathbf{0}. Since any basis for a subspace of \mathbb{R}^n extends to a basis of the whole space, choose vectors \mathbf{v}_{k+1},\dots,\mathbf{v}_n so that \{\mathbf{v}_1,\dots,\mathbf{v}_n\} is a basis for \mathbb{R}^n. The claim is that the n-k image vectors A\mathbf{v}_{k+1},\dots,A\mathbf{v}_n form a basis for \operatorname{Col}A, which pins \operatorname{rank}(A) = n-k and finishes the proof.

Spanning. Every \mathbf{b}\in\operatorname{Col}A has the form A\mathbf{x} for some \mathbf{x}\in\mathbb{R}^n. Write \mathbf{x} in the extended basis as \mathbf{x} = c_1\mathbf{v}_1 + \cdots + c_n\mathbf{v}_n. Applying A and using linearity, \mathbf{b} = A\mathbf{x} = c_1A\mathbf{v}_1 + \cdots + c_kA\mathbf{v}_k + c_{k+1}A\mathbf{v}_{k+1} + \cdots + c_nA\mathbf{v}_n. The first k terms vanish because \mathbf{v}_1,\dots,\mathbf{v}_k \in \operatorname{Nul}A, leaving \mathbf{b} = c_{k+1}A\mathbf{v}_{k+1} + \cdots + c_nA\mathbf{v}_n. So the n-k image vectors already span \operatorname{Col}A.

Independence. Suppose c_{k+1}A\mathbf{v}_{k+1} + \cdots + c_nA\mathbf{v}_n = \mathbf{0}. By linearity this is A(c_{k+1}\mathbf{v}_{k+1} + \cdots + c_n\mathbf{v}_n) = \mathbf{0}, which puts the vector \mathbf{w} = c_{k+1}\mathbf{v}_{k+1} + \cdots + c_n\mathbf{v}_n inside \operatorname{Nul}A. Since \{\mathbf{v}_1,\dots,\mathbf{v}_k\} is a basis for \operatorname{Nul}A, \mathbf{w} also equals some combination d_1\mathbf{v}_1 + \cdots + d_k\mathbf{v}_k. Subtracting the two expressions for \mathbf{w} gives d_1\mathbf{v}_1 + \cdots + d_k\mathbf{v}_k - c_{k+1}\mathbf{v}_{k+1} - \cdots - c_n\mathbf{v}_n = \mathbf{0}, a linear relation among all n vectors of the basis \{\mathbf{v}_1,\dots,\mathbf{v}_n\}. Basis vectors are linearly independent, so every coefficient here is 0 — in particular c_{k+1} = \cdots = c_n = 0. Only the trivial combination of A\mathbf{v}_{k+1},\dots,A\mathbf{v}_n gives \mathbf{0}, so these vectors are linearly independent.

Together, spanning and independence make \{A\mathbf{v}_{k+1},\dots, A\mathbf{v}_n\} a basis for \operatorname{Col}A, so \operatorname{rank}(A) = \dim\left(\operatorname{Col}A\right) = n-k. Moving k to the other side recovers \operatorname{rank}(A) + \dim\left(\operatorname{Nul}A\right) = n.

Sanity check: for the 3\times 3 identity matrix, \operatorname{Nul}A = \{\mathbf{0}\}, so k = 0 and \operatorname{rank}(A) = 3 = n, matching the fact that the identity is invertible and its columns already span \mathbb{R}^3. At the other extreme, the 3\times 3 zero matrix has \operatorname{Nul}A = \mathbb{R}^3, so k=3 and \operatorname{rank}(A) = 0 — no columns to span anything — and 0 + 3 = 3 still checks out.

Start from what makes T linear: T(\mathbf{u}+\mathbf{v}) = T(\mathbf{u}) + T(\mathbf{v}) and T(c\mathbf{u}) = cT(\mathbf{u}) for every scalar c and every \mathbf{u},\mathbf{v}\in\mathbb{R}^n. Repeated use of both rules extends linearity from two vectors to any finite combination: T(c_1\mathbf{u}_1 + c_2\mathbf{u}_2 + \cdots + c_k\mathbf{u}_k) = c_1T(\mathbf{u}_1) + c_2T(\mathbf{u}_2) + \cdots + c_kT(\mathbf{u}_k).

Now take any \mathbf{x}\in\mathbb{R}^n and write it in coordinates, \mathbf{x} = (x_1,x_2,\dots,x_n). The standard basis vectors \mathbf{e}_1,\dots,\mathbf{e}_n (each \mathbf{e}_i having a 1 in slot i and 0 elsewhere) reconstruct \mathbf{x} exactly as the combination \mathbf{x} = x_1\mathbf{e}_1 + x_2\mathbf{e}_2 + \cdots + x_n\mathbf{e}_n. Apply T to both sides and invoke the extended linearity rule above: T(\mathbf{x}) = x_1T(\mathbf{e}_1) + x_2T(\mathbf{e}_2) + \cdots + x_nT(\mathbf{e}_n).

This right-hand side is a weighted sum of the fixed vectors T(\mathbf{e}_1),\dots,T(\mathbf{e}_n), each in \mathbb{R}^m, weighted by the entries of \mathbf{x}. That is precisely what matrix-vector multiplication computes when those vectors are placed as the columns of a matrix A: by the column definition of the matrix-vector product, A\mathbf{x} = x_1A_{:,1} + x_2A_{:,2} + \cdots + x_nA_{:,n}, where A_{:,i} denotes column i of A. Matching this term by term against the expansion of T(\mathbf{x}) above forces A_{:,i} = T(\mathbf{e}_i) for every i — so setting A = \begin{bmatrix} T(\mathbf{e}_1) & T(\mathbf{e}_2) & \cdots & T(\mathbf{e}_n) \end{bmatrix} gives T(\mathbf{x}) = A\mathbf{x} for every \mathbf{x}\in\mathbb{R}^n, as claimed. A is m\times n because each column T(\mathbf{e}_i) lives in \mathbb{R}^m and there are n of them, one per input coordinate.

Uniqueness: suppose some other matrix B also satisfies T(\mathbf{x}) = B\mathbf{x} for every \mathbf{x}. Taking \mathbf{x} = \mathbf{e}_i gives B\mathbf{e}_i = T(\mathbf{e}_i) — but B\mathbf{e}_i is just column i of B, by the same column definition of matrix-vector multiplication used above. So column i of B equals T(\mathbf{e}_i) for every i, which means B = A. No other matrix represents T.

Sanity check: take T:\mathbb{R}^2\to\mathbb{R}^2 the 90° counterclockwise rotation, so T(\mathbf{e}_1) = (0,1) and T(\mathbf{e}_2) = (-1,0). The theorem says the standard matrix is A = \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix}. Testing it on \mathbf{x} = (3,1): A\mathbf{x} = (-1,3), which is exactly (3,1) rotated 90° counterclockwise — matching T directly, as it must.

To show \ker T is a subspace of V, three things must hold: it contains the zero vector, it is closed under addition, and it is closed under scalar multiplication.

Contains the zero vector. Linearity forces T(\mathbf{0}) = \mathbf{0}: writing \mathbf{0} = 0\mathbf{v} for any \mathbf{v} \in V, T(\mathbf{0}) = T(0\mathbf{v}) = 0\, T(\mathbf{v}) = \mathbf{0}. So \mathbf{0} \in \ker T, and \ker T is nonempty.

Closed under addition. Let \mathbf{u}, \mathbf{v} \in \ker T, so T(\mathbf{u}) = \mathbf{0} and T(\mathbf{v}) = \mathbf{0}. Then T(\mathbf{u} + \mathbf{v}) = T(\mathbf{u}) + T(\mathbf{v}) = \mathbf{0} + \mathbf{0} = \mathbf{0}, using additivity of T. Hence \mathbf{u} + \mathbf{v} \in \ker T.

Closed under scalar multiplication. Let \mathbf{v} \in \ker T and let c be any scalar. Then T(c\mathbf{v}) = c\, T(\mathbf{v}) = c\mathbf{0} = \mathbf{0}, using homogeneity of T. Hence c\mathbf{v} \in \ker T.

All three closure properties hold, so \ker T is a subspace of V.

Sanity check: if T is the zero transformation, T(\mathbf{v}) = \mathbf{0} for every \mathbf{v} \in V, so \ker T = V itself — the whole space, which is certainly a subspace of V. At the other extreme, if T is one-to-one, the only vector sent to \mathbf{0} is \mathbf{0} itself, so \ker T = \{\mathbf{0}\}, the trivial subspace. Both extremes are genuine subspaces, exactly as the theorem guarantees.

(\Rightarrow) Suppose T is one-to-one. Linearity already guarantees T(\mathbf{0}) = \mathbf{0}, so \mathbf{0} \in \ker(T) no matter what. The question is whether anything else can join it. Suppose some vector \mathbf{v} \in \ker(T), so T(\mathbf{v}) = \mathbf{0}. Since we also know T(\mathbf{0}) = \mathbf{0}, this vector \mathbf{v} and the zero vector map to the same output. One-to-one means distinct inputs never share an output, so \mathbf{v} = \mathbf{0}. No vector other than \mathbf{0} can survive in \ker(T), which forces \ker(T) = \{\mathbf{0}\}.

(\Leftarrow) Suppose \ker(T) = \{\mathbf{0}\}. To show T is one-to-one, take any two vectors \mathbf{u}, \mathbf{v} \in V with T(\mathbf{u}) = T(\mathbf{v}), and show \mathbf{u} = \mathbf{v}. Subtract: T(\mathbf{u}) - T(\mathbf{v}) = \mathbf{0}. Linearity lets the left side collapse into a single application of T to a difference of vectors: T(\mathbf{u} - \mathbf{v}) = \mathbf{0}. This says \mathbf{u} - \mathbf{v} \in \ker(T). But the kernel holds only the zero vector, so \mathbf{u} - \mathbf{v} = \mathbf{0} \quad\Longrightarrow\quad \mathbf{u} = \mathbf{v}. Equal outputs forced equal inputs, which is exactly the definition of one-to-one.

Why this trick works at all: for a general function, checking "one-to-one" means comparing every pair of inputs to see if they collide. For a linear transformation, that whole pairwise search collapses onto a single set — the kernel — because any collision T(\mathbf{u}) = T(\mathbf{v}) can always be rewritten as T(\mathbf{u} - \mathbf{v}) = \mathbf{0}, turning a statement about two vectors into a statement about one. That's the move this proof makes in both directions, and it's the reason injectivity of a linear map is a one-vector test instead of an all-pairs test.

Since B is a basis, every v \in V has a unique representation v = c_1b_1 + \cdots + c_nb_n, so the coordinate vector [v]_B = (c_1, \dots, c_n) is well defined and T is a well-defined function from V to \mathbb{R}^n. Three things remain: T is linear, injective, and surjective.

Linearity. Take u, v \in V with [u]_B = (a_1,\dots,a_n) and [v]_B = (c_1,\dots,c_n), and take a scalar k. Then u + v = \sum_{i=1}^n a_ib_i + \sum_{i=1}^n c_ib_i = \sum_{i=1}^n (a_i+c_i)b_i, so [u+v]_B = (a_1+c_1,\dots,a_n+c_n) = [u]_B + [v]_B, i.e. T(u+v) = T(u) + T(v). Likewise kv = \sum_{i=1}^n (ka_i)b_i gives [kv]_B = k[v]_B, i.e. T(kv) = kT(v). So T is linear.

Injective. Suppose T(v) = 0, i.e. [v]_B = (0,\dots,0). Then v = 0b_1 + \cdots + 0b_n = 0, so \ker T = \{0\}. A linear map with trivial kernel is injective: if T(u) = T(v) then linearity gives T(u-v) = 0, forcing u - v = 0, i.e. u = v.

Surjective. Let (c_1,\dots,c_n) \in \mathbb{R}^n be arbitrary, and set v = c_1b_1 + \cdots + c_nb_n \in V. Since the b_i form a basis, this sum is the unique representation of v in terms of B, so T(v) = [v]_B = (c_1,\dots,c_n). Every vector in \mathbb{R}^n is hit.

T is a linear bijection, so it is an isomorphism, and V \cong \mathbb{R}^n.

Sanity check: take V = P_2 (polynomials of degree at most 2) with basis B = \{1, x, x^2\}, and \dim V = 3. The polynomial 5 - 2x + 7x^2 maps to (5,-2,7) \in \mathbb{R}^3, addition and scaling of polynomials match componentwise addition and scaling of their coefficient triples, and every triple in \mathbb{R}^3 comes from exactly one polynomial in P_2 — exactly the linear, injective, surjective behavior the theorem predicts.

Every \mathbf{x} \in V has a unique expansion in the \mathcal{B}-basis, \mathbf{x} = x_1\mathbf{b}_1 + x_2\mathbf{b}_2 + \cdots + x_n\mathbf{b}_n, where the scalars x_1, \dots, x_n are exactly the entries of [\mathbf{x}]_{\mathcal B}. Take \mathcal{C}-coordinates of both sides. Coordinate-taking relative to a fixed basis is linear — the coordinate map sends sums to sums and scalar multiples to scalar multiples, since it just reads off the unique coefficients in a linear combination — so [\mathbf{x}]_{\mathcal C} = [x_1\mathbf{b}_1 + \cdots + x_n\mathbf{b}_n]_{\mathcal C} = x_1[\mathbf{b}_1]_{\mathcal C} + x_2[\mathbf{b}_2]_{\mathcal C} + \cdots + x_n[\mathbf{b}_n]_{\mathcal C}. The right side is precisely a linear combination of the columns of P_{\mathcal C \leftarrow \mathcal B} with weights x_1, \dots, x_n, which is the definition of the matrix-vector product: x_1[\mathbf{b}_1]_{\mathcal C} + \cdots + x_n[\mathbf{b}_n]_{\mathcal C} = \begin{bmatrix} [\mathbf{b}_1]_{\mathcal C} & \cdots & [\mathbf{b}_n]_{\mathcal C} \end{bmatrix}\begin{bmatrix} x_1 \\ \vdots \\ x_n \end{bmatrix} = P_{\mathcal C \leftarrow \mathcal B}\,[\mathbf{x}]_{\mathcal B}. Combining the two displays gives [\mathbf{x}]_{\mathcal C} = P_{\mathcal C \leftarrow \mathcal B}\,[\mathbf{x}]_{\mathcal B}, as claimed.

For invertibility, run the same argument with the roles of \mathcal{B} and \mathcal{C} swapped: the matrix P_{\mathcal B \leftarrow \mathcal C} built from the \mathcal{B}-coordinates of the \mathcal{C}-basis vectors satisfies [\mathbf{x}]_{\mathcal B} = P_{\mathcal B \leftarrow \mathcal C}\,[\mathbf{x}]_{\mathcal C} for every \mathbf{x} \in V. Substituting the first formula into the second (and vice versa) gives [\mathbf{x}]_{\mathcal B} = P_{\mathcal B \leftarrow \mathcal C}P_{\mathcal C \leftarrow \mathcal B}\,[\mathbf{x}]_{\mathcal B}, \qquad [\mathbf{x}]_{\mathcal C} = P_{\mathcal C \leftarrow \mathcal B}P_{\mathcal B \leftarrow \mathcal C}\,[\mathbf{x}]_{\mathcal C} for every \mathbf{x} \in V. A matrix that fixes every coordinate vector in \mathbb{R}^n must be the identity matrix I_n (test it against the standard basis vectors, whose coordinate images run over all of \mathbb{R}^n), so both products equal I_n. Hence P_{\mathcal C \leftarrow \mathcal B} and P_{\mathcal B \leftarrow \mathcal C} are inverse to one another.

Columns are orthonormal. Write Q = \begin{bmatrix} q_1 & q_2 & \cdots & q_n \end{bmatrix} in terms of its columns. The (i,j) entry of Q^{\mathsf T}Q is the dot product q_i \cdot q_j, so the defining equation Q^{\mathsf T}Q = I_n says exactly q_i \cdot q_j = \begin{cases} 1 & i = j \\ 0 & i \neq j \end{cases}, i.e. each column has unit length and distinct columns are orthogonal — an orthonormal set. Running the argument with QQ^{\mathsf T} in place of Q^{\mathsf T}Q shows the rows are orthonormal too, once the inverse identity below is in hand.

Inverse equals transpose. Q^{\mathsf T}Q = I_n is already the statement that Q^{\mathsf T} is a left inverse of Q. Since Q is square, a left inverse is a two-sided inverse: Q^{\mathsf T}Q = I_n forces Q to have full rank n, so Q is invertible, and multiplying Q^{\mathsf T}Q = I_n on the right by Q^{-1} gives Q^{\mathsf T} = Q^{-1}. In particular QQ^{\mathsf T} = QQ^{-1} = I_n as well.

Determinant is \pm 1. Take determinants of both sides of Q^{\mathsf T}Q = I_n: \det(Q^{\mathsf T})\det(Q) = \det(I_n) = 1. A matrix and its transpose share the same determinant, so this reads (\det Q)^2 = 1, and therefore \det Q = \pm 1.

Length is preserved. For any vector x, \lVert Qx\rVert^2 = (Qx)\cdot(Qx) = (Qx)^{\mathsf T}(Qx) = x^{\mathsf T}Q^{\mathsf T}Qx = x^{\mathsf T}I_nx = x^{\mathsf T}x = \lVert x\rVert^2, using Q^{\mathsf T}Q = I_n in the middle step. Taking square roots (both sides are nonnegative) gives \lVert Qx\rVert = \lVert x\rVert: multiplying by Q never stretches or shrinks a vector, so Q acts as a rotation, reflection, or some composition of the two.

Write \hat{\mathbf v} = \operatorname{proj}_W(\mathbf v) for the orthogonal projection of \mathbf v onto W, and recall the defining property of that projection: the residual \mathbf v - \hat{\mathbf v} is orthogonal to every vector in W, i.e. \mathbf v - \hat{\mathbf v} \in W^\perp.

Let \mathbf w be any vector in W with \mathbf w \ne \hat{\mathbf v}. Split the error \mathbf v - \mathbf w into two pieces by inserting and subtracting \hat{\mathbf v}: \mathbf v - \mathbf w = (\mathbf v - \hat{\mathbf v}) + (\hat{\mathbf v} - \mathbf w). The first piece, \mathbf v - \hat{\mathbf v}, lies in W^\perp by the projection property above. The second piece, \hat{\mathbf v} - \mathbf w, is a difference of two vectors already in W, so it lies in W itself (subspaces are closed under subtraction). These two pieces are therefore orthogonal to each other, and the Pythagorean theorem for orthogonal vectors gives \lVert\mathbf v - \mathbf w\rVert^2 = \lVert\mathbf v - \hat{\mathbf v}\rVert^2 + \lVert\hat{\mathbf v} - \mathbf w\rVert^2.

Since \mathbf w \ne \hat{\mathbf v}, the second term \lVert\hat{\mathbf v} - \mathbf w\rVert^2 is strictly positive, so \lVert\mathbf v - \mathbf w\rVert^2 > \lVert\mathbf v - \hat{\mathbf v}\rVert^2, and taking square roots (both sides are nonnegative) gives \lVert\mathbf v - \mathbf w\rVert > \lVert\mathbf v - \hat{\mathbf v}\rVert. Because \mathbf w was an arbitrary point of W other than \hat{\mathbf v}, no other point of W can match or beat the distance achieved by \hat{\mathbf v} — so \hat{\mathbf v} is the unique closest point of W to \mathbf v.

Existence. Let \{\mathbf{u}_1, \dots, \mathbf{u}_p\} be an orthogonal basis for W (one always exists, by Gram–Schmidt if no other is in hand). Define

Each term of \hat{\mathbf y} is a scalar multiple of some \mathbf{u}_i \in W, and W is closed under linear combinations, so \hat{\mathbf y} \in W outright — no further check needed there. The work is showing \mathbf{z} \in W^\perp, and for that it suffices to show \mathbf{z} is orthogonal to every basis vector \mathbf{u}_i, since orthogonality to a spanning set of W forces orthogonality to all of W. Fix i and compute \mathbf{z}\cdot\mathbf{u}_i, using that \mathbf{u}_i\cdot\mathbf{u}_j = 0 for j \ne i (orthogonality of the basis) so every term of \hat{\mathbf y} except the i-th drops out of the dot product:

So \mathbf{z} \perp \mathbf{u}_i for every i, hence \mathbf{z} \in W^\perp, and \mathbf{y} = \hat{\mathbf y} + \mathbf{z} is a decomposition of the required form. That \hat{\mathbf y} is exactly \operatorname{proj}_W\mathbf{y} is just the definition of the projection onto W built from an orthogonal basis of W.

Uniqueness. Suppose \mathbf{y} also decomposes as \mathbf{y} = \hat{\mathbf y}' + \mathbf{z}' with \hat{\mathbf y}' \in W and \mathbf{z}' \in W^\perp. Subtracting the two decompositions of \mathbf{y} gives

The left side lies in W, since W is a subspace and both \hat{\mathbf y}, \hat{\mathbf y}' \in W; the right side lies in W^\perp for the same reason applied to \mathbf{z}, \mathbf{z}' \in W^\perp. So this common vector belongs to W \cap W^\perp. But a vector orthogonal to every vector in W — including itself, if it lies in W — must have zero length: taking \mathbf{v} = \hat{\mathbf y} - \hat{\mathbf y}' as that vector, \mathbf{v} \in W and \mathbf{v} \in W^\perp together force \mathbf{v}\cdot\mathbf{v} = 0, so \mathbf{v} = \mathbf{0}. Hence \hat{\mathbf y} = \hat{\mathbf y}' and, from the subtracted equation, \mathbf{z} = \mathbf{z}' as well. The decomposition is unique.

Sanity check: take W = \operatorname{span}\{(1,0,0),(0,1,0)\}, the xy-plane inside \mathbb{R}^3, and \mathbf{y} = (2,3,5). The orthogonal basis \{(1,0,0),(0,1,0)\} gives \hat{\mathbf y} = (2,3,0) and \mathbf{z} = (0,0,5), which is indeed orthogonal to both basis vectors — exactly the geometrically obvious split of \mathbf{y} into its in-plane part and its part along the perpendicular z-axis.

The claim has two parts: each \mathbf{u}_k built this way is orthogonal to every \mathbf{u}_i that came before it, and the running set \{\mathbf{u}_1, \ldots, \mathbf{u}_k\} spans the same subspace as \{\mathbf{v}_1, \ldots, \mathbf{v}_k\} at every stage. Both are proved together by induction on k.

Base case. \mathbf{u}_1 = \mathbf{v}_1, so trivially \operatorname{span}\{\mathbf{u}_1\} = \operatorname{span}\{\mathbf{v}_1\}, and there is no earlier vector to be orthogonal to.

Inductive step. Suppose \{\mathbf{u}_1, \ldots, \mathbf{u}_{k-1}\} is already orthogonal and spans the same set as \{\mathbf{v}_1, \ldots, \mathbf{v}_{k-1}\}. Form \mathbf{u}_k = \mathbf{v}_k - \sum_{i=1}^{k-1} \operatorname{proj}_{\mathbf{u}_i}(\mathbf{v}_k) = \mathbf{v}_k - \sum_{i=1}^{k-1} \frac{\mathbf{v}_k \cdot \mathbf{u}_i}{\mathbf{u}_i \cdot \mathbf{u}_i}\,\mathbf{u}_i.

Orthogonality. Fix any j with 1 \le j \le k-1 and dot both sides with \mathbf{u}_j: \mathbf{u}_k \cdot \mathbf{u}_j = \mathbf{v}_k \cdot \mathbf{u}_j - \sum_{i=1}^{k-1} \frac{\mathbf{v}_k \cdot \mathbf{u}_i}{\mathbf{u}_i \cdot \mathbf{u}_i}\,(\mathbf{u}_i \cdot \mathbf{u}_j). By the inductive hypothesis \{\mathbf{u}_1, \ldots, \mathbf{u}_{k-1}\} is orthogonal, so \mathbf{u}_i \cdot \mathbf{u}_j = 0 for every i \ne j in that sum, and every term except i = j drops out. What's left is \mathbf{u}_k \cdot \mathbf{u}_j = \mathbf{v}_k \cdot \mathbf{u}_j - \frac{\mathbf{v}_k \cdot \mathbf{u}_j}{\mathbf{u}_j \cdot \mathbf{u}_j}\,(\mathbf{u}_j \cdot \mathbf{u}_j) = \mathbf{v}_k \cdot \mathbf{u}_j - \mathbf{v}_k \cdot \mathbf{u}_j = 0. Since j was arbitrary, \mathbf{u}_k is orthogonal to every earlier \mathbf{u}_i, and since \{\mathbf{u}_1, \ldots, \mathbf{u}_{k-1}\} is already pairwise orthogonal, \{\mathbf{u}_1, \ldots, \mathbf{u}_k\} is too. (This step needs \mathbf{u}_j \ne \mathbf{0} for each earlier j, which the spanning argument below guarantees: a zero \mathbf{u}_j would mean \mathbf{v}_j already lay in the span of the vectors before it, contradicting linear independence of the original basis.)

Same span. The formula for \mathbf{u}_k writes \mathbf{v}_k as \mathbf{u}_k plus a linear combination of \mathbf{u}_1, \ldots, \mathbf{u}_{k-1}, i.e. \mathbf{v}_k = \mathbf{u}_k + \sum_{i=1}^{k-1} \frac{\mathbf{v}_k \cdot \mathbf{u}_i}{\mathbf{u}_i \cdot \mathbf{u}_i}\,\mathbf{u}_i. By the inductive hypothesis each \mathbf{u}_i on the right is itself a linear combination of \mathbf{v}_1, \ldots, \mathbf{v}_{k-1}, so \mathbf{v}_k \in \operatorname{span}\{\mathbf{v}_1, \ldots, \mathbf{v}_{k-1}, \mathbf{u}_k\}, and conversely \mathbf{u}_k \in \operatorname{span}\{\mathbf{v}_1, \ldots, \mathbf{v}_k\} directly from its defining formula. Combined with the inductive hypothesis that \operatorname{span}\{\mathbf{u}_1, \ldots, \mathbf{u}_{k-1}\} = \operatorname{span}\{\mathbf{v}_1, \ldots, \mathbf{v}_{k-1}\}, this gives \operatorname{span}\{\mathbf{u}_1, \ldots, \mathbf{u}_k\} = \operatorname{span}\{\mathbf{v}_1, \ldots, \mathbf{v}_k\}.

By induction, both properties hold for every k up to n, so \{\mathbf{u}_1, \ldots, \mathbf{u}_n\} is an orthogonal basis for the same subspace spanned by \{\mathbf{v}_1, \ldots, \mathbf{v}_n\}. None of the \mathbf{u}_k can be \mathbf{0}: that would force \mathbf{v}_k into the span of \mathbf{v}_1, \ldots, \mathbf{v}_{k-1}, contradicting linear independence, so every division by \mathbf{u}_i \cdot \mathbf{u}_i above is division by a positive number.

Normalizing. Once \{\mathbf{u}_1, \ldots, \mathbf{u}_n\} is orthogonal with no zero vector among them, setting \mathbf{q}_k = \mathbf{u}_k / \lVert\mathbf{u}_k\rVert scales each vector to unit length without disturbing the directions, so \{\mathbf{q}_1, \ldots, \mathbf{q}_n\} is orthonormal and spans the same subspace.

Run the Gram–Schmidt process on the columns \mathbf{a}_1, \dots, \mathbf{a}_n, producing orthonormal vectors \mathbf{q}_1, \dots, \mathbf{q}_n one at a time. At each step, build an intermediate vector \mathbf{w}_k by stripping off \mathbf{a}_k's component along every \mathbf{q}_i already found, \mathbf{w}_k = \mathbf{a}_k - \sum_{i=1}^{k-1} (\mathbf{a}_k \cdot \mathbf{q}_i)\,\mathbf{q}_i, \qquad \mathbf{q}_k = \frac{\mathbf{w}_k}{\lVert\mathbf{w}_k\rVert}. Linear independence of \mathbf{a}_1, \dots, \mathbf{a}_n guarantees \mathbf{w}_k \neq \mathbf{0} at every step (a zero \mathbf{w}_k would make \mathbf{a}_k a combination of \mathbf{a}_1, \dots, \mathbf{a}_{k-1}), so each division is legal and \lVert\mathbf{w}_k\rVert > 0.

Solve the defining relation for \mathbf{a}_k itself by moving the sum to the other side and substituting \mathbf{w}_k = \lVert\mathbf{w}_k\rVert\, \mathbf{q}_k: \mathbf{a}_k = \sum_{i=1}^{k-1} (\mathbf{a}_k \cdot \mathbf{q}_i)\, \mathbf{q}_i \;+\; \lVert\mathbf{w}_k\rVert\,\mathbf{q}_k. Every vector on the right is one of \mathbf{q}_1, \dots, \mathbf{q}_k, so this says each column of A is a linear combination of the \mathbf{q}_i's with an explicit, computable set of coefficients — exactly the data needed to assemble a matrix equation.

Let Q be the m \times n matrix with columns \mathbf{q}_1, \dots, \mathbf{q}_n, and define R entry by entry from the coefficients just found: r_{ik} = \mathbf{a}_k \cdot \mathbf{q}_i \ \ (i < k), \qquad r_{kk} = \lVert\mathbf{w}_k\rVert, \qquad r_{ik} = 0 \ \ (i > k). The zero pattern is exactly upper triangular, since \mathbf{a}_k was only ever built from \mathbf{q}_1, \dots, \mathbf{q}_k — it carries no component along \mathbf{q}_{k+1}, \dots, \mathbf{q}_n, which is what r_{ik} = 0 for i > k records. With this R, column k of the product QR is \sum_{i=1}^{n} r_{ik}\,\mathbf{q}_i = \sum_{i=1}^{k} r_{ik}\,\mathbf{q}_i, which is precisely the expression for \mathbf{a}_k derived above. Since this holds for every column, A = QR.

Finally, Q has orthonormal columns by construction of Gram–Schmidt, so Q^{\mathsf T}Q = I_n, and R is upper triangular with diagonal entries r_{kk} = \lVert\mathbf{w}_k\rVert > 0, as required.

Sanity check: take n = 1, so A is a single nonzero column \mathbf{a}_1. Then \mathbf{w}_1 = \mathbf{a}_1, \mathbf{q}_1 = \mathbf{a}_1 / \lVert\mathbf{a}_1\rVert, and R = [\lVert\mathbf{a}_1\rVert] — a single positive number. Indeed QR = \mathbf{q}_1 \cdot \lVert\mathbf{a}_1\rVert = \mathbf{a}_1, recovering A exactly, as the factorization must in the smallest nontrivial case.

A least-squares solution \hat{\mathbf x} is, by definition, a minimizer of \lVert A\mathbf x - \mathbf b \rVert over all \mathbf x — equivalently, a vector for which A\hat{\mathbf x} is the point of \operatorname{Col}(A) closest to \mathbf b. The closest-point characterization of a projection says exactly that the residual \mathbf b - A\hat{\mathbf x} must be orthogonal to \operatorname{Col}(A): no vector in that column space can be nearer to \mathbf b than its orthogonal projection, so if the residual had any component lying in \operatorname{Col}(A), subtracting that component would move A\hat{\mathbf x} closer to \mathbf b still, contradicting minimality.

Orthogonality to every vector in \operatorname{Col}(A) means orthogonality to every column of A in particular, since those columns span \operatorname{Col}(A). Writing \mathbf a_1, \dots, \mathbf a_n for the columns of A, this is the system of dot-product conditions

Stacking these n scalar equations into one vector equation is precisely what A^{\mathsf T} does: the j-th row of A^{\mathsf T} is \mathbf a_j^{\mathsf T}, so left-multiplying by A^{\mathsf T} produces a vector whose j-th entry is \mathbf a_j \cdot (\mathbf b - A\hat{\mathbf x}). Since every one of those entries is 0,

Distributing A^{\mathsf T} across the difference and moving the A^{\mathsf T} A \hat{\mathbf x} term to the other side gives

Sanity check: if \mathbf b already lies in \operatorname{Col}(A), the system A\mathbf x = \mathbf b is consistent and any exact solution \mathbf x^* leaves residual \mathbf b - A\mathbf x^* = \mathbf 0, which is trivially orthogonal to every column of A — so \mathbf x^* satisfies the normal equations too, and least squares reduces to ordinary solving exactly when there's no error left to minimize.

Write \mathbf{x} for a candidate vector and split the target \mathbf{b} against the subspace \operatorname{Col}(A): every \mathbf{b} decomposes uniquely as \mathbf{b} = \hat{\mathbf{b}} + \mathbf{b}^{\perp}, where \hat{\mathbf{b}} = \operatorname{proj}_{\operatorname{Col}(A)}\mathbf{b} is the closest point to \mathbf{b} inside \operatorname{Col}(A) and \mathbf{b}^{\perp} is orthogonal to every vector in \operatorname{Col}(A).

For any \mathbf{x}, the vector A\mathbf{x} lies in \operatorname{Col}(A), so the residual splits as A\mathbf{x} - \mathbf{b} = (A\mathbf{x} - \hat{\mathbf{b}}) - \mathbf{b}^{\perp}, a sum of a vector inside \operatorname{Col}(A) and a vector orthogonal to it. The Pythagorean theorem for orthogonal vectors then gives \lVert A\mathbf{x} - \mathbf{b} \rVert^2 = \lVert A\mathbf{x} - \hat{\mathbf{b}} \rVert^2 + \lVert \mathbf{b}^{\perp} \rVert^2. The second term does not depend on \mathbf{x} at all, so this quantity is smallest exactly when the first term is 0 — that is, exactly when A\mathbf{x} = \hat{\mathbf{b}}. Every minimizer \hat{\mathbf{x}} therefore satisfies A\hat{\mathbf{x}} = \hat{\mathbf{b}}, which is equivalent to \mathbf{b} - A\hat{\mathbf{x}} = \mathbf{b} - \hat{\mathbf{b}} = \mathbf{b}^{\perp} \in \operatorname{Col}(A)^{\perp}.

A vector lies in \operatorname{Col}(A)^{\perp} exactly when it is orthogonal to every column of A, i.e. when A^{\mathsf T} sends it to \mathbf{0} — because each entry of A^{\mathsf T}\mathbf{v} is one such dot product. So the residual condition above rewrites as A^{\mathsf T}\left(\mathbf{b} - A\hat{\mathbf{x}}\right) = \mathbf{0} \quad\Longleftrightarrow\quad A^{\mathsf T} A \hat{\mathbf{x}} = A^{\mathsf T} \mathbf{b}, the normal equations. Every step above is reversible, so a vector \hat{\mathbf{x}} satisfies the normal equations if and only if it minimizes \lVert A\mathbf{x} - \mathbf{b} \rVert.

For uniqueness under linear independence: if the columns of A are linearly independent then A\mathbf{v} = \mathbf{0} forces \mathbf{v} = \mathbf{0}, so for any nonzero \mathbf{v}, \mathbf{v}^{\mathsf T} A^{\mathsf T} A \mathbf{v} = (A\mathbf{v})^{\mathsf T}(A\mathbf{v}) = \lVert A\mathbf{v} \rVert^2 > 0, which shows A^{\mathsf T} A \mathbf{v} \neq \mathbf{0} for every nonzero \mathbf{v} — so A^{\mathsf T} A has trivial null space and, being square, is invertible. Multiplying both sides of the normal equations by \left(A^{\mathsf T} A\right)^{-1} then isolates the unique solution \hat{\mathbf{x}} = \left(A^{\mathsf T} A\right)^{-1} A^{\mathsf T} \mathbf{b}.

Run the Gram–Schmidt process on the columns \mathbf{a}_1, \dots, \mathbf{a}_n, producing orthonormal vectors \mathbf{q}_1, \dots, \mathbf{q}_n one at a time. At each step, build an intermediate vector \mathbf{w}_k by stripping off \mathbf{a}_k's component along every \mathbf{q}_i already found, \mathbf{w}_k = \mathbf{a}_k - \sum_{i=1}^{k-1} (\mathbf{a}_k \cdot \mathbf{q}_i)\,\mathbf{q}_i, \qquad \mathbf{q}_k = \frac{\mathbf{w}_k}{\lVert\mathbf{w}_k\rVert}. Linear independence of \mathbf{a}_1, \dots, \mathbf{a}_n guarantees \mathbf{w}_k \neq \mathbf{0} at every step (a zero \mathbf{w}_k would make \mathbf{a}_k a combination of \mathbf{a}_1, \dots, \mathbf{a}_{k-1}), so each division is legal and \lVert\mathbf{w}_k\rVert > 0.

Solve the defining relation for \mathbf{a}_k itself by moving the sum to the other side and substituting \mathbf{w}_k = \lVert\mathbf{w}_k\rVert\, \mathbf{q}_k: \mathbf{a}_k = \sum_{i=1}^{k-1} (\mathbf{a}_k \cdot \mathbf{q}_i)\, \mathbf{q}_i \;+\; \lVert\mathbf{w}_k\rVert\,\mathbf{q}_k. Every vector on the right is one of \mathbf{q}_1, \dots, \mathbf{q}_k, so this says each column of A is a linear combination of the \mathbf{q}_i's with an explicit, computable set of coefficients — exactly the data needed to assemble a matrix equation.

Let Q be the m \times n matrix with columns \mathbf{q}_1, \dots, \mathbf{q}_n, and define R entry by entry from the coefficients just found: r_{ik} = \mathbf{a}_k \cdot \mathbf{q}_i \ \ (i < k), \qquad r_{kk} = \lVert\mathbf{w}_k\rVert, \qquad r_{ik} = 0 \ \ (i > k). The zero pattern is exactly upper triangular, since \mathbf{a}_k was only ever built from \mathbf{q}_1, \dots, \mathbf{q}_k — it carries no component along \mathbf{q}_{k+1}, \dots, \mathbf{q}_n, which is what r_{ik} = 0 for i > k records. With this R, column k of the product QR is \sum_{i=1}^{n} r_{ik}\,\mathbf{q}_i = \sum_{i=1}^{k} r_{ik}\,\mathbf{q}_i, which is precisely the expression for \mathbf{a}_k derived above. Since this holds for every column, A = QR.

Finally, Q has orthonormal columns by construction of Gram–Schmidt, so Q^{\mathsf T}Q = I_n, and R is upper triangular with diagonal entries r_{kk} = \lVert\mathbf{w}_k\rVert > 0, as required.

Sanity check: take n = 1, so A is a single nonzero column \mathbf{a}_1. Then \mathbf{w}_1 = \mathbf{a}_1, \mathbf{q}_1 = \mathbf{a}_1 / \lVert\mathbf{a}_1\rVert, and R = [\lVert\mathbf{a}_1\rVert] — a single positive number. Indeed QR = \mathbf{q}_1 \cdot \lVert\mathbf{a}_1\rVert = \mathbf{a}_1, recovering A exactly, as the factorization must in the smallest nontrivial case.

Since the columns of A are linearly independent, A^{\mathsf T}A is invertible and the least-squares solution is the unique solution of the normal equations A^{\mathsf T}A\,\mathbf{x} = A^{\mathsf T}\mathbf{b}.

Substitute A = QR into both sides. On the left, A^{\mathsf T}A = (QR)^{\mathsf T}(QR) = R^{\mathsf T}Q^{\mathsf T}QR. The columns of Q are orthonormal, which is exactly the statement Q^{\mathsf T}Q = I_n, so this collapses to A^{\mathsf T}A = R^{\mathsf T}R. On the right, A^{\mathsf T}\mathbf{b} = (QR)^{\mathsf T}\mathbf{b} = R^{\mathsf T}Q^{\mathsf T}\mathbf{b}. The normal equations become R^{\mathsf T}R\,\mathbf{x} = R^{\mathsf T}Q^{\mathsf T}\mathbf{b}.

R is invertible — it is a factor of the full-column-rank matrix A in a factorization where Q has full column rank — so R^{\mathsf T} is also invertible, and left-multiplying both sides by \left(R^{\mathsf T}\right)^{-1} cancels it entirely: R\mathbf{x} = Q^{\mathsf T}\mathbf{b}. Since R is invertible this system has the unique solution \hat{\mathbf{x}} = R^{-1}Q^{\mathsf T}\mathbf{b}, which must therefore be the least-squares solution guaranteed by the normal equations.

Why this is the practical version of the normal equations. Forming A^{\mathsf T}A squares the conditioning of the problem — round-off error that was small in A gets amplified when A is multiplied against itself. The QR route never forms A^{\mathsf T}A at all: it computes Q^{\mathsf T}\mathbf{b} directly and then solves the triangular system R\hat{\mathbf x} = Q^{\mathsf T}\mathbf{b} by back substitution, which is both cheaper and numerically safer. The two formulas agree exactly — R^{\mathsf T}R above is nothing but a triangular factorization of A^{\mathsf T}A itself — they just reach the same \hat{\mathbf{x}} by different arithmetic paths.

By definition, \lambda is an eigenvalue of A exactly when there exists a nonzero vector \mathbf{x} with A\mathbf{x} = \lambda\mathbf{x}. Move everything to one side. Since \lambda\mathbf{x} = \lambda I\mathbf{x}, this rearranges to A\mathbf{x} - \lambda I\mathbf{x} = \mathbf{0}, \qquad\text{i.e.}\qquad (A - \lambda I)\mathbf{x} = \mathbf{0}. So \lambda is an eigenvalue of A precisely when the homogeneous system (A - \lambda I)\mathbf{x} = \mathbf{0} has a solution \mathbf{x} other than the trivial one — equivalently, precisely when \text{Nul}(A - \lambda I) contains more than just \mathbf{0}.

A square matrix has a nontrivial null space if and only if it fails to be invertible: if A - \lambda I were invertible, multiplying both sides of (A - \lambda I)\mathbf{x} = \mathbf{0} by (A - \lambda I)^{-1} would force \mathbf{x} = \mathbf{0}, ruling out a nonzero solution. So a nonzero \mathbf{x} can exist only when A - \lambda I is singular.

Finally, a square matrix is singular exactly when its determinant is zero. Chaining the equivalences together: \lambda \text{ is an eigenvalue of } A \iff (A-\lambda I)\mathbf{x} = \mathbf{0} \text{ has a nontrivial solution} \iff \det(A - \lambda I) = 0.

Sanity check: take A = I, the identity matrix. Then A - \lambda I = (1-\lambda)I, and \det\left((1-\lambda)I\right) = (1-\lambda)^n, which vanishes only at \lambda = 1. That matches the identity matrix's only eigenvalue: every nonzero vector satisfies I\mathbf{x} = 1\cdot\mathbf{x}, so \lambda = 1 is the sole eigenvalue, exactly as the characteristic equation predicts.

The eigenvalues of A are the scalars \lambda for which A - \lambda I is singular, i.e. the roots of the characteristic polynomial p(\lambda) = \det(A - \lambda I). Subtracting \lambda I from A only touches the diagonal — every off-diagonal entry of A - \lambda I equals the corresponding off-diagonal entry of A — so A - \lambda I is triangular of the same shape as A, with diagonal entries a_{11} - \lambda, a_{22} - \lambda, \dots, a_{nn} - \lambda.

Take A upper triangular; the lower-triangular and diagonal cases run the same way. Expand \det(A - \lambda I) by cofactors down the first column. Because every entry below the (1,1) position is 0, the only surviving cofactor is the one for the (1,1) entry: \det(A - \lambda I) = (a_{11} - \lambda)\,\det\!\left(M_{11}\right), where M_{11} is the (n-1)\times(n-1) matrix left after deleting row 1 and column 1 of A - \lambda I. Deleting the first row and column of an upper triangular matrix leaves an upper triangular matrix, and its diagonal is a_{22} - \lambda, \dots, a_{nn} - \lambda — exactly the diagonal of A - \lambda I with the first entry removed.

Repeating the same first-column expansion on M_{11} peels off the factor (a_{22} - \lambda) and leaves an (n-2)\times(n-2) upper triangular matrix with diagonal a_{33} - \lambda, \dots, a_{nn} - \lambda. Continuing down to a single 1 \times 1 block (a_{nn} - \lambda) gives p(\lambda) = \det(A - \lambda I) = (a_{11} - \lambda)(a_{22} - \lambda) \cdots (a_{nn} - \lambda).

This is already a product of linear factors, so its roots are read off directly: p(\lambda) = 0 exactly when \lambda = a_{ii} for some i, and a diagonal entry that repeats k times among a_{11}, \dots, a_{nn} contributes a factor (a_{ii} - \lambda)^k, giving that eigenvalue algebraic multiplicity k. Hence the eigenvalues of A are precisely its diagonal entries, with multiplicity.

Sanity check: take A = I_n, the identity matrix — diagonal, with every diagonal entry equal to 1. The theorem says every eigenvalue of I_n is 1, matching \det(I_n - \lambda I) = (1-\lambda)^n, which is correct since I_n \mathbf{x} = \mathbf{x} = 1 \cdot \mathbf{x} for every nonzero \mathbf{x}.

The eigenvalues of A are, by definition, the roots of the characteristic polynomial p(\lambda) = \det(\lambda I - A), counted with the same multiplicity they have as roots. Since p is a monic polynomial of degree n (this comes out of the Leibniz expansion below), it factors completely over \mathbb{C} as

Multiplying this out and collecting by degree gives

where the middle terms are omitted because only the leading coefficient, the \lambda^{n-1} coefficient, and the constant term will be needed.

The constant term recovers \det(A). Setting \lambda = 0 in the factored form gives p(0) = \prod_{i=1}^n(-\lambda_i) = (-1)^n\prod_{i=1}^n \lambda_i, matching the constant term already written in (1). On the other hand, p(0) = \det(0\cdot I - A) = \det(-A), and pulling a factor of -1 out of each of the n rows of -A gives \det(-A) = (-1)^n\det(A). So (1)'s constant term is also (-1)^n\det(A). Comparing the two expressions for p(0), (-1)^n\det(A) = (-1)^n\prod_{i=1}^n \lambda_i \quad\Longrightarrow\quad \det(A) = \prod_{i=1}^n \lambda_i.

The \lambda^{n-1} coefficient recovers \operatorname{tr}(A). Write out p(\lambda) = \det(\lambda I - A) by the Leibniz formula, summing over every permutation \sigma of \{1,\dots,n\}: p(\lambda) = \sum_{\sigma} \operatorname{sgn}(\sigma)\prod_{i=1}^n (\lambda I - A)_{i,\sigma(i)}. An off-diagonal entry of \lambda I - A, where \sigma(i)\ne i, is just -a_{i,\sigma(i)} — it carries no \lambda at all. Only a diagonal entry, where \sigma(i) = i, contributes the factor (\lambda - a_{ii}), which is linear in \lambda. So the degree in \lambda of the term for a given \sigma equals the number of fixed points of \sigma. Any permutation other than the identity moves at least two indices, so it fixes at most n-2 points and contributes a term of degree at most n-2 in \lambda — too low to affect the \lambda^n or \lambda^{n-1} coefficients. Those two top coefficients therefore come entirely from the identity permutation's term, \prod_{i=1}^n(\lambda - a_{ii}) = \lambda^n - \left(\sum_{i=1}^n a_{ii}\right)\lambda^{n-1} + \cdots = \lambda^n - \operatorname{tr}(A)\,\lambda^{n-1} + \cdots, using \operatorname{tr}(A) = \sum_i a_{ii}. (This also confirms p is monic of degree n, as claimed above.) Comparing this \lambda^{n-1} coefficient, -\operatorname{tr}(A), with the one already recorded in (1), -\sum_i \lambda_i, gives -\operatorname{tr}(A) = -\sum_{i=1}^n \lambda_i \quad\Longrightarrow\quad \operatorname{tr}(A) = \sum_{i=1}^n \lambda_i.

Sanity check: for a 2\times 2 matrix A = \begin{pmatrix} a & b \\ c & d\end{pmatrix} this says \lambda_1+\lambda_2 = a+d and \lambda_1\lambda_2 = ad-bc — exactly the sum and product of roots of the familiar quadratic \lambda^2 - (a+d)\lambda + (ad-bc) = 0 obtained directly from \det(\lambda I - A) = 0, which is the n=2 case of the argument above.

(\Leftarrow) n independent eigenvectors give A = PDP^{-1}. Suppose v_1, \dots, v_n are linearly independent eigenvectors of A, with Av_i = \lambda_i v_i. Let P be the matrix with these vectors as its columns, and let D = \operatorname{diag}(\lambda_1, \dots, \lambda_n). Multiplying A into P column by column, AP = A\begin{bmatrix} v_1 & \cdots & v_n \end{bmatrix} = \begin{bmatrix} Av_1 & \cdots & Av_n \end{bmatrix} = \begin{bmatrix} \lambda_1 v_1 & \cdots & \lambda_n v_n \end{bmatrix}. Multiplying P by D on the right scales each column of P by the matching diagonal entry of D, producing the very same matrix: PD = \begin{bmatrix} v_1 & \cdots & v_n \end{bmatrix} \operatorname{diag}(\lambda_1,\dots,\lambda_n) = \begin{bmatrix} \lambda_1 v_1 & \cdots & \lambda_n v_n \end{bmatrix}. So AP = PD. Because P's columns are n linearly independent vectors in \mathbb{R}^n, P is invertible, and right-multiplying both sides by P^{-1} gives A = PDP^{-1}.

(\Rightarrow) A = PDP^{-1} gives n independent eigenvectors. Suppose A = PDP^{-1} with P invertible and D = \operatorname{diag}(\lambda_1, \dots, \lambda_n). Then AP = PD, and reading that equation column by column exactly as above shows Av_i = \lambda_i v_i, where v_i is the i-th column of P. Each v_i is nonzero — a zero column would make P singular — so v_1,\dots,v_n are eigenvectors of A. Since P is invertible, its columns are linearly independent, which produces the needed n independent eigenvectors.

Sanity check: take A = \begin{bmatrix} 4 & 1 \\ 0 & 3 \end{bmatrix}, upper triangular, so its eigenvalues sit on the diagonal: \lambda_1 = 4, \lambda_2 = 3. Solving (A - 4I)v = 0 gives eigenvector v_1 = (1,0), and (A-3I)v = 0 gives v_2 = (1,-1). These two are not scalar multiples of each other, so they're independent, and the theorem guarantees A = PDP^{-1} with P = \begin{bmatrix} 1 & 1 \\ 0 & -1 \end{bmatrix} and D = \begin{bmatrix} 4 & 0 \\ 0 & 3 \end{bmatrix} — which a direct multiplication PDP^{-1} confirms reproduces A.

Characteristic polynomial (and hence eigenvalues). Write the characteristic polynomial of B using B = P^{-1}AP: \det(B - \lambda I) = \det\!\left(P^{-1}AP - \lambda I\right) = \det\!\left(P^{-1}AP - \lambda P^{-1}P\right) = \det\!\left(P^{-1}(A - \lambda I)P\right). The determinant is multiplicative, so this splits into three factors: \det(B - \lambda I) = \det(P^{-1})\,\det(A - \lambda I)\,\det(P) = \det(A - \lambda I), since \det(P^{-1}) = 1/\det(P) cancels against \det(P). So A and B share the same characteristic polynomial \det(A - \lambda I), and therefore the same roots — the same eigenvalues, each with the same algebraic multiplicity, since multiplicity is read off the same polynomial.

Determinant. Setting \lambda = 0 in the identity just proved gives \det(B) = \det(A) directly. (Equivalently, \det(B) = \det(P^{-1})\det(A)\det(P) = \det(A) by the same multiplicative argument.)

Trace. The characteristic polynomial of an n \times n matrix M has the form \det(M - \lambda I) = (-\lambda)^n + \operatorname{tr}(M)(-\lambda)^{n-1} + \cdots + \det(M), so \operatorname{tr}(M) is exactly the coefficient of (-\lambda)^{n-1}. Since A and B share the same characteristic polynomial, they share that coefficient, so \operatorname{tr}(A) = \operatorname{tr}(B).

Rank. Left-multiplication by the invertible matrix P^{-1} and right-multiplication by the invertible matrix P each preserve rank — multiplying by an invertible matrix on either side sends a basis of the column space to another basis of a column space of the same dimension, and never collapses two independent columns together, since invertible maps are one-to-one. So \operatorname{rank}(B) = \operatorname{rank}(P^{-1}AP) = \operatorname{rank}(AP) = \operatorname{rank}(A).

Why this matters. These quantities are exactly the ones that don't depend on which basis is used to write a linear transformation as a matrix — changing basis replaces A by P^{-1}AP for the change-of-basis matrix P, so \det, \operatorname{tr}, \operatorname{rank}, and the eigenvalues are properties of the underlying transformation itself, not artifacts of coordinates. The converse is false: matrices can share all of these invariants without being similar (a diagonalizable matrix is never similar to a non-diagonalizable one with the same eigenvalues), so this theorem gives necessary conditions for similarity, not sufficient ones.

Argue by contradiction, and suppose the set is linearly dependent. Among all dependence relations c_1 v_1 + c_2 v_2 + \cdots + c_k v_k = \mathbf{0} with not every c_i zero, pick one using the fewest nonzero vectors — say it uses exactly m of them, relabeled so they're v_1, \ldots, v_m with every coefficient c_1, \ldots, c_m nonzero: c_1 v_1 + c_2 v_2 + \cdots + c_m v_m = \mathbf{0}. \qquad (\star) Since eigenvectors are nonzero by definition, m \geq 2 — a single nonzero vector can't satisfy c_1 v_1 = \mathbf{0} with c_1 \neq 0.

Apply A to both sides of (\star) and use Av_i = \lambda_i v_i: c_1 \lambda_1 v_1 + c_2 \lambda_2 v_2 + \cdots + c_m \lambda_m v_m = \mathbf{0}. \qquad (\star\star) Now multiply (\star) by \lambda_m and subtract (\star\star) from it, so the v_m terms cancel: c_1(\lambda_m - \lambda_1) v_1 + c_2(\lambda_m - \lambda_2) v_2 + \cdots + c_{m-1}(\lambda_m - \lambda_{m-1}) v_{m-1} = \mathbf{0}. This is a dependence relation on only m - 1 of the vectors — fewer than the minimal count m assumed above. For it to be genuinely shorter (not just relabeled with a zero coefficient tacked on), at least one of its coefficients must be nonzero. Since the eigenvalues are distinct, \lambda_m - \lambda_i \neq 0 for every i < m, and c_i \neq 0 by choice of (\star), so each coefficient c_i(\lambda_m - \lambda_i) is nonzero for i = 1, \ldots, m-1. That contradicts m being the fewest vectors any dependence relation among v_1, \ldots, v_k can use.

The contradiction shows no dependence relation exists, so \{v_1, \ldots, v_k\} is linearly independent.

Sanity check: take k = 2. The argument reduces to: if c_1 v_1 + c_2 v_2 = \mathbf{0} with c_1 \neq 0, applying A gives c_1\lambda_1 v_1 + c_2\lambda_2 v_2 = \mathbf{0}; subtracting \lambda_2 times the first equation leaves c_1(\lambda_1 - \lambda_2)v_1 = \mathbf{0}, and since \lambda_1 \neq \lambda_2 and v_1 \neq \mathbf{0}, this forces c_1 = 0 — contradicting c_1 \neq 0. So two eigenvectors for distinct eigenvalues are always independent, exactly the base case the induction-style argument above builds on.

(\Leftarrow) Orthogonal diagonalizability forces symmetry. Suppose A = QDQ^{\mathsf T} with Q orthogonal and D diagonal. Diagonal matrices equal their own transpose, and transposing a product reverses the order of its factors, so A^{\mathsf T} = (QDQ^{\mathsf T})^{\mathsf T} = (Q^{\mathsf T})^{\mathsf T}D^{\mathsf T}Q^{\mathsf T} = QDQ^{\mathsf T} = A. So A is symmetric. The rest of the proof handles the harder direction.

Symmetric matrices have only real eigenvalues. Let \lambda be an eigenvalue of A, allowing for the moment that it might be complex, with eigenvector v (possibly with complex entries), so Av = \lambda v. Let \bar v denote the entrywise complex conjugate of v, and consider the scalar \bar v^{\mathsf T} A v. Computing it one way, \bar v^{\mathsf T} A v = \bar v^{\mathsf T}(\lambda v) = \lambda (\bar v^{\mathsf T} v) = \lambda \lVert v \rVert^2, since \bar v^{\mathsf T} v = \sum_i \overline{v_i} v_i = \sum_i |v_i|^2 = \lVert v \rVert^2 is real and positive (v \neq 0, being an eigenvector). Computing the same scalar the other way, using A^{\mathsf T} = A and that \bar v^{\mathsf T} A v is a 1\times 1 matrix and so equal to its own transpose, \bar v^{\mathsf T} A v = \left(\bar v^{\mathsf T} A v\right)^{\mathsf T} = v^{\mathsf T} A^{\mathsf T} \bar v = v^{\mathsf T} A \bar v. Taking the complex conjugate of Av = \lambda v gives A\bar v = \bar\lambda \bar v (since A has real entries, conjugation passes through it untouched), so v^{\mathsf T} A \bar v = v^{\mathsf T}(\bar\lambda \bar v) = \bar\lambda (v^{\mathsf T}\bar v) = \bar\lambda \lVert v\rVert^2, using v^{\mathsf T}\bar v = \overline{\bar v^{\mathsf T} v} = \lVert v\rVert^2 again. The two computations of \bar v^{\mathsf T} A v must agree: \lambda \lVert v \rVert^2 = \bar\lambda \lVert v \rVert^2. Since \lVert v \rVert^2 \neq 0, this forces \lambda = \bar\lambda, i.e. \lambda is real. So every eigenvalue of A is real, and consequently every eigenvalue has a real eigenvector (a real matrix minus a real multiple of I_n has a real null space).

Eigenvectors for distinct eigenvalues are orthogonal. Let Av_1 = \lambda_1 v_1 and Av_2 = \lambda_2 v_2 with \lambda_1 \neq \lambda_2. Using symmetry, \lambda_1 (v_1 \cdot v_2) = (Av_1) \cdot v_2 = v_1^{\mathsf T}A^{\mathsf T} v_2 = v_1^{\mathsf T}Av_2 = v_1 \cdot (Av_2) = \lambda_2(v_1 \cdot v_2). So (\lambda_1 - \lambda_2)(v_1 \cdot v_2) = 0, and since \lambda_1 \neq \lambda_2, this forces v_1 \cdot v_2 = 0.

Building Q by induction on n. The base case n=1 is trivial: any 1\times 1 matrix is already diagonal, and Q = [1] is orthogonal. For the inductive step, assume every symmetric (n-1)\times(n-1) matrix is orthogonally diagonalizable, and let A be a symmetric n \times n matrix. By the real-eigenvalue fact above, A has a real eigenvalue \lambda_1 with a real, nonzero eigenvector; rescale it to a unit vector q_1, so Aq_1 = \lambda_1 q_1 and \lVert q_1 \rVert = 1. Extend \{q_1\} to an orthonormal basis \{q_1, q_2, \dots, q_n\} of \mathbb{R}^n (extend to any basis and apply Gram–Schmidt, which leaves q_1 fixed since it is already a unit vector orthogonal to nothing else yet chosen), and let P be the orthogonal matrix with these vectors as columns.

The subspace W = \operatorname{span}(q_2,\dots,q_n) = q_1^\perp is invariant under A: for w \in W, (Aw)\cdot q_1 = w \cdot (Aq_1) = w \cdot (\lambda_1 q_1) = \lambda_1(w\cdot q_1) = 0, using symmetry in the first step exactly as above, so Aw is again orthogonal to q_1, i.e. Aw \in W. Consequently P^{\mathsf T}AP has the block form P^{\mathsf T}AP = \begin{bmatrix} \lambda_1 & 0 \\ 0 & A_1 \end{bmatrix} for some symmetric (n-1)\times(n-1) matrix A_1: the first column of P^{\mathsf T}AP is P^{\mathsf T}(Aq_1) = P^{\mathsf T}(\lambda_1 q_1) = \lambda_1 e_1 (since P^{\mathsf T}q_1 = e_1, q_1 being the first column of the orthogonal matrix P), giving the zero block below \lambda_1; and P^{\mathsf T}AP is symmetric because A is, which forces the zero block to the right of \lambda_1 as well. By the inductive hypothesis, A_1 = Q_1D_1Q_1^{\mathsf T} for some orthogonal (n-1)\times(n-1) matrix Q_1 and diagonal D_1. Setting Q = P\begin{bmatrix} 1 & 0 \\ 0 & Q_1 \end{bmatrix}, \qquad D = \begin{bmatrix} \lambda_1 & 0 \\ 0 & D_1 \end{bmatrix}, the block-diagonal factor is orthogonal (its columns are orthonormal, being e_1 together with the orthonormal columns of Q_1), so Q is a product of orthogonal matrices and hence orthogonal itself. Unwinding the block product, Q^{\mathsf T}AQ = \begin{bmatrix} 1 & 0 \\ 0 & Q_1^{\mathsf T} \end{bmatrix}P^{\mathsf T}AP\begin{bmatrix} 1 & 0 \\ 0 & Q_1 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & Q_1^{\mathsf T} \end{bmatrix}\begin{bmatrix} \lambda_1 & 0 \\ 0 & A_1 \end{bmatrix}\begin{bmatrix} 1 & 0 \\ 0 & Q_1 \end{bmatrix} = \begin{bmatrix} \lambda_1 & 0 \\ 0 & Q_1^{\mathsf T}A_1Q_1 \end{bmatrix} = D, so A = QDQ^{\mathsf T} with Q orthogonal and D diagonal, completing the induction.

Sanity check: take the symmetric matrix A = \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}. Its characteristic equation is (2-\lambda)^2 - 1 = 0, so \lambda_1 = 3 and \lambda_2 = 1 — both real, as guaranteed. Solving (A-3I)v=0 gives v_1 = (1,1), and (A-I)v=0 gives v_2 = (1,-1); these two are orthogonal (v_1 \cdot v_2 = 1 - 1 = 0), exactly as the theorem predicts for eigenvectors of distinct eigenvalues. Normalizing gives q_1 = \frac{1}{\sqrt2}(1,1) and q_2 = \frac{1}{\sqrt2}(1,-1), and setting Q = \begin{bmatrix} 1/\sqrt2 & 1/\sqrt2 \\ 1/\sqrt2 & -1/\sqrt2 \end{bmatrix}, D = \begin{bmatrix} 3 & 0 \\ 0 & 1 \end{bmatrix}, a direct multiplication confirms QDQ^{\mathsf T} = A.

Step 1: \lambda = 1 is an eigenvalue of P. Let \mathbf{1}^{\mathsf T} be the row vector of all ones. Each column of P sums to 1, and summing a column is exactly the dot product of that column with \mathbf{1}^{\mathsf T}, so \mathbf{1}^{\mathsf T} P = \mathbf{1}^{\mathsf T}. This says \mathbf{1} is a left eigenvector of P for eigenvalue 1, i.e. \mathbf{1}^{\mathsf T}(P - I) = \mathbf{0}^{\mathsf T}, so (P-I)^{\mathsf T} has a nontrivial null vector and is therefore singular. A matrix and its transpose share the same determinant, so \det(P - I) = \det\big((P-I)^{\mathsf T}\big) = 0 as well — meaning P - I is singular in its own right, which makes 1 an eigenvalue of P itself (with an honest eigenvector \mathbf{q}, not just a left one).

Step 2: every other eigenvalue is strictly smaller in modulus. This is the step that uses regularity. Because some power P^k has strictly positive entries, no proper subset of coordinates can be permuted among itself forever under repeated multiplication by P — every coordinate eventually mixes with every other one. That mixing rules out a second eigenvalue on the unit circle: two independent eigenvectors sharing |\lambda| = 1 would force P^k to keep some coordinate pattern rigid no matter how large k gets, contradicting the fact that P^k's columns are eventually all strictly positive and hence free to blend. So among the eigenvalues of P, \lambda = 1 is the unique one on the unit circle, and every other eigenvalue \lambda_i satisfies |\lambda_i| < 1.

Step 3: convergence. Suppose for the moment that P is diagonalizable, with eigenvalues 1 = \lambda_1, \lambda_2, \dots, \lambda_n (repeated by multiplicity) and corresponding eigenvectors \mathbf{v}_1 = \mathbf{q}, \mathbf{v}_2, \dots, \mathbf{v}_n forming a basis of \mathbb{R}^n. Write any starting probability vector in that basis: \mathbf{x}_0 = c_1\mathbf{q} + c_2\mathbf{v}_2 + \cdots + c_n\mathbf{v}_n. Applying P repeatedly scales each eigencomponent by its own eigenvalue: P^k\mathbf{x}_0 = c_1\mathbf{q} + c_2\lambda_2^k\mathbf{v}_2 + \cdots + c_n\lambda_n^k\mathbf{v}_n. By Step 2, |\lambda_i| < 1 for every i \ge 2, so \lambda_i^k \to 0 as k \to \infty, and every term but the first vanishes in the limit: P^k\mathbf{x}_0 \to c_1\mathbf{q}. To pin down c_1, apply \mathbf{1}^{\mathsf T} to both sides of the eigenbasis expansion. Multiplying a probability vector by P never changes its total: \mathbf{1}^{\mathsf T}P\mathbf{x} = \mathbf{1}^{\mathsf T}\mathbf{x} by Step 1's identity, so \mathbf{1}^{\mathsf T}P^k\mathbf{x}_0 = \mathbf{1}^{\mathsf T}\mathbf{x}_0 = 1 for every k. Taking k\to\infty on the left gives \mathbf{1}^{\mathsf T}(c_1\mathbf{q}) = 1, and since \mathbf{q} is normalized to have entries summing to 1, this forces c_1 = 1. Hence P^k\mathbf{x}_0 \to \mathbf{q} for every probability vector \mathbf{x}_0 — the limit is the same vector \mathbf{q} no matter where the chain started, which is exactly what makes \mathbf{q} the steady state rather than one of many. (When P is not diagonalizable the same conclusion holds via its Jordan form instead of an eigenbasis: the \lambda = 1 block contributes the constant term c_1\mathbf{q} exactly as above, while every other block carries a factor \lambda_i^k or k^m\lambda_i^k that still decays to 0 since |\lambda_i| < 1.)

Uniqueness of \mathbf{q} among steady-state vectors now follows for free: if \mathbf{q}' were another vector with P\mathbf{q}' = \mathbf{q}' and entries summing to 1, then \mathbf{q}' = P^k\mathbf{q}' \to \mathbf{q} for every k, so \mathbf{q}' = \mathbf{q}. Positivity of \mathbf{q} follows the same way, taking \mathbf{x}_0 to be any single standard basis vector \mathbf{e}_j: P^k\mathbf{e}_j is a column of P^k, which is strictly positive once k is large enough, and its limit \mathbf{q} inherits that positivity (a limit of positive vectors can only fail to be positive by having some entry equal to 0 in the limit, but each entry of \mathbf{q} is itself the limit of a sequence bounded below by the smallest entry of the eventually-positive P^k, which stays bounded away from 0).