We are going to establish the an identity for the Vandermonde determinants:

\left|\matrix{
1 & a_{1} & a_{1}^{2} & \ldots & a_{1}^{n-1} \\
1 & a_{2} & a_{2}^{2} & \ldots & a_{2}^{n-1} \\
1 & a_{3} & a_{3}^{2} & \ldots & a_{3}^{n-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & a_{n} & a_{n}^{2} & \ldots & a_{n}^{n-1} \\
}\right| = \prod_{1\leq i<j\leq n}(a_{j}-a_{i})
.

The proof is by induction on the size n of the determinant. For n=1 the claim is true by default. The first "real" example comes with the next value n=2:

\left|\matrix{
1 & a_{1} \\
1 & a_{2} \\
}\right| = a_{2}-a_{1}
.

For the inductive step, assume

\left|\matrix{
1 & a_{1} & a_{1}^{2} & \ldots & a_{1}^{k-1} \\
1 & a_{2} & a_{2}^{2} & \ldots & a_{2}^{k-1} \\
1 & a_{3} & a_{3}^{2} & \ldots & a_{3}^{k-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & a_{k} & a_{k}^{2} & \ldots & a_{k}^{k-1} \\
}\right| = \prod_{1\leq i<j\leq k}(a_{j}-a_{i})
.

and let n=k+1. Consider

\Delta=\left|\matrix{
1 & a_{1} & a_{1}^{2} & \ldots & a_{1}^{k} \\
1 & a_{2} & a_{2}^{2} & \ldots & a_{2}^{k} \\
1 & a_{3} & a_{3}^{2} & \ldots & a_{3}^{k} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & a_{k+1} & a_{k+1}^{2} & \ldots & a_{k+1}^{k} \\
}\right|
.

First subtract the last row in turn from all others:

\Delta=\left|\matrix{
0 & a_{1}-a_{k+1} & a_{1}^{2}-a_{k+1}^{2} & \ldots & a_{1}^{k-1}-a_{k+1}^{k-1} & a_{1}^{k}-a_{k+1}^{k} \\
0 & a_{2}-a_{k+1} & a_{2}^{2}-a_{k+1}^{2} & \ldots & a_{2}^{k-1}-a_{k+1}^{k-1} & a_{2}^{k}-a_{k+1}^{k} \\
0 & a_{3}-a_{k+1} & a_{3}^{2}-a_{k+1}^{2} & \ldots & a_{3}^{k-1}-a_{k+1}^{k-1} & a_{3}^{k}-a_{k+1}^{k} \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
1 & a_{k+1} & a_{k+1}^{2} & \ldots & a_{k+1}^{k-1} & a_{k+1}^{k} \\
}\right|
.

Next, multiply column k by a_{k+1} and subtract the result from column k+1; multiply column k-1 by a_{k} and subtract the result from column k, and so on:

\Delta=\left|\matrix{
0 & a_{1}-a_{k+1} & a_{1}(a_{1}-a_{k+1}) & \ldots & a_{1}^{k-2}(a_{1}-a_{k+1}) & a_{1}^{k-1}(a_{1}-a_{k+1}) \\
0 & a_{2}-a_{k+1} & a_{2}(a_{2}-a_{k+1}) & \ldots & a_{2}^{k-2}(a_{2}-a_{k+1}) & a_{2}^{k-1}(a_{2}-a_{k+1}) \\
0 & a_{3}-a_{k+1} & a_{3}(a_{3}-a_{k+1}) & \ldots & a_{3}^{k-2}(a_{3}-a_{k+1}) & a_{3}^{k-1}(a_{3}-a_{k+1}) \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
1 & 0 & 0 & \ldots & 0 & 0 \\
}\right|
.

Expanding the determinant by the first column gives

\Delta=(-1)^{k+2}\prod_{1\leq i\leq k}(a_{i}-a_{k+1})\left|\matrix{
1 & a_{1} & a_{1}^{2} & \ldots & a_{1}^{k-1} \\
1 & a_{2} & a_{2}^{2} & \ldots & a_{2}^{k-1} \\
1 & a_{3} & a_{3}^{2} & \ldots & a_{3}^{k-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & a_{k} & a_{k}^{2} & \ldots & a_{k}^{k-1} \\
}\right|=
\prod_{1\leq i\leq k}(a_{k+1}-a_{i})\left|\matrix{
1 & a_{1} & a_{1}^{2} & \ldots & a_{1}^{k-1} \\
1 & a_{2} & a_{2}^{2} & \ldots & a_{2}^{k-1} \\
1 & a_{3} & a_{3}^{2} & \ldots & a_{3}^{k-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & a_{k} & a_{k}^{2} & \ldots & a_{k}^{k-1} \\
}\right|
,

which - together with the inductive assumption - completes the proof.

Alternate Solution from Properties of Determinants

Since the determinant |M| of a matrix M is a sum of products of entries chosen one from each column and row, it is a polynomial in the variables a_1, \cdots, a_n, and its degree is 1+2+\cdots+n-1 = {n\choose 2}. Since a_i = a_j makes |M|=0, the polynomial must have roots such that it is divisible by (a_j-a_i), and there are {n\choose 2} such factors. Since this matches the degree of the polynomial, we must have all of its factors; hence the polynomial must be C\prod_{1\leq i<j\leq n}(a_{j}-a_{i}) for some constant C.

To determine C, note that one term in the determinant is the product of the diagonal elements, which occurs with coefficient +1 and has the property that a_n has the highest possible power, a_{n-1} the next highest, and so on. In expanding the polynomial product, the powers can only be maximized in this way if we prefer, in each factor, the term with the largest index. Since the form of the sum shows that this is always the positive term, we find that the coefficient of the diagonal product is C = +1, and the result follows.