One fine summer morning not so long ago, I woke up with a start to the stark realization that I had no idea how to derive the formula for the sum of powers of natural numbers. I am talking about a closed-form expression for the series (for non-negative integer \(p\)) \[ \sum_{i=1}^n {i^p} = 1 + {2^p} + {3^p} + \dots + {n^p} \]
Assuming that we know that the formula is a polynomial of degree \((p+1)\) (which happens to be the case), we can easily derive the coefficients by substituting for \(n=1,2,\dots p+1\), and solving the resulting simultaneous linear equations. But what if we do not know that?
Rather than look for an answer online immediately, I decided to kill some time thinking about the problem by myself. I was able to reach a solution without too much effort. As a cursory search did not reveal the exact method described elsewhere, I am sharing my derivation here.
Let us denote the sum \(\sum_{i=1}^n i^p\) as \(S_p(n)\).
Case \(p=0\)
We will take our base case as \(p=0\). Obviously, \(\sum_{i=1}^n i^0 = 1 + 1 + \dots + 1\) is just \(n\).
Case \(p=1\)
To compute the formula for \(p=1\), consider the following sum:
\begin{align} T_0 &= (1^0) + (1^0+2^0) + \dots + (1^0 + 2^0 + \dots + n^0) \end{align}
Now, \(T_0\) is the same as \(1+2+\dots + n = S_1(n)\). But as there are \(n\) number of \(1^0\), \((n-1)\) number of \(2^0\), etc., we can write \(T_0\) also in the following manner:
\begin{align}
T_0 &=1^0 . n + 2^0 . (n-1) + \dots n^0 . (n - (n-1))\\\
&= \sum_{i=1}^{n} (i^0 . (n - i + 1))\\\
&=\sum_{i=1}^{n} ((n+1) - i)\\\
&=\sum_{i=1}^{n} (i^0(n+1) - i)\\\
&=(n+1)\sum_{i=1}^n i^0 - \sum_{i=1}^{n} i\\\
&=(n+1)S_0(n) - S_1(n)
\end{align}
Now we can equate the two expressions for \(T_0\), and get:
\begin{align}
S_1(n) &= (n+1)S_0(n) - S_1(n)\\\
\therefore 2S_1(n) &= (n+1)S_0(n)\\\
\therefore S_1(n) &= \frac{1}{2}(n+1)S_0(n)\\\
&= \frac{1}{2}n(n+1)
\end{align}
Which checks out!
Case \(p=2\)
\begin{align}
T_1 &= (1^1) + (1^1+2^1) + \dots + (1^1 + 2^1 + \dots + n^1) \\\
&=\sum_{i=1}^n S_1(i) \\\
&=\sum_{i=1}^n \frac{1}{2}i(i+1)\\\
&=\sum_{i=1}^n \frac{1}{2}(i^2+i)\\\
&=\frac{1}{2}\sum_{i=1}^n i^2 + \frac{1}{2}\sum_{i=1}^n i\\\
&=\frac{1}{2} S_2(n) + \frac{1}{2} S_1(n)
\end{align}
Rewriting \(T_1\) as in the previous case, we see that
\begin{align}
T_1 &=1^1 . n + 2^1 . (n-1) + \dots n^1 . (n - (n-1))\\\
&= \sum_{i=1}^{n} (i^1 . (n - i + 1))\\\
&=\sum_{i=1}^{n} (i(n+1) - i^2)\\\
&=\sum_{i=1}^{n} (i^1(n+1) - i^2)\\\
&=(n+1)\sum_{i=1}^n i^1 - \sum_{i=1}^{n} i^2\\\
&=(n+1)S_1(n) - S_2(n)
\end{align}
Equating, we can obtain the formula for \(S_2(n)\).
\begin{align}
\frac{1}{2} S_2(n) + \frac{1}{2} S_1(n) &= (n+1)S_1(n) - S_2(n) \\\
\therefore \frac{3}{2} S_2(n) &= (n+\frac{1}{2}) S_1(n) \\\
\therefore S_2(n) &= \frac{2}{3} (n+\frac{1}{2}) S_1(n) \\\
&= \frac{1}{3} (2n + 1) \frac{n (n+1)}{2}\\\
&= \frac{1}{6} n(n+1)(2n+1)
\end{align}
General Case
Now we can generalize what we have been doing so far. But first, we need to assume that \(S_p(n)\) is a polynomial of degree \((p+1)\) with no constant term. But wait! Isn’t that cheating? In fact it is not, because we are introducing this assumption as part of an induction hypothesis. In other words, assuming that \(S_p(n)\) is a polynomial of degree \((p+1)\), we will prove that \(S_{p+1}(n)\) is similarly a polynomial of degree \((p+2)\) with no constant term. We will also get the complete (recursive) formula for the sum of powers of natural numbers.
So let us assume that \(S_p(n) = a_{p,1}n^{p+1} + a_{p,2}n^{p} + \dots + a_{p, p+1}n\), for all \(p \le k\). We can take one of \(p=0, 1, 2\) as the base case.
For computing the formula for \(p=k+1\), consider the expression:
\begin{align}
T_{k} &= (1^k) + (1^k+2^k) + \dots + (1^k + 2^k + \dots + n^k)\\\
&= \sum_{i=1}^{n} S_k(i)\\\
&= a_{k,1} S_{k+1}(n) + a_{k,2} S_k(n) + \dots + a_{k, k+1} S_1(n)
\end{align}
But also,
\begin{align}
T_k &= 1^k . n + 2^k . (n-1) + \dots + n^k . (n - (n-1))\\\
&= \sum_{i=1}^{n} i^k . (n-i+1)\\\
&= (n+1) \sum_{i=1}^{n} i^k - \sum_{i=1}^{n} i^{k+1}\\\
&= (n+1) S_{k}(n) - S_{k+1}(n)
\end{align}
Equating both the expressions, we get
\begin{align}
&a_{k,1} S_{k+1}(n) + a_{k,2} S_k(n) + \dots + a_{k, k+1} S_1(n) &\\\
&,,,=(n +1)
S_{k}(n) -
S_{k+1}(n)\\\
&(a_{k,1} + 1)S_{k+1}(n) \\\
&,,,= (n+1 - a_{k,2})S_k(n) - a_{k,3} S_{k-1}(n) - \dots
- a_{k,k+1} S_1(n)\\\
&,,,= (n+1 - a_{k,2})S_k(n) - \sum_{j=3}^{k+1} a_{k,j}S_{k+2-j}(n)\\\
&,,,= nS_k(n) - (a_{k,2} - 1)S_k(n) - \sum_{j=3}^{k+1} a_{k,j}S_{k+2-j}(n)
\end{align}
From the last expression we can see that the largest power of \(S_{k+1}(n)\) is \(n\) times the largest power of \(S_k(n)\). Since \(S_k(n)\) is a polynomial of degree \(k+1\), \(S_{k+1}(n)\) must be a polynomial of degree \((k+2)\). Also note that all the coefficients of the recursive sum are just the coefficients of \(S_k(n)\). So if we have \(S_1, S_2, \dots S_k(n)\), we can readily compute the expression for \(S_{k+1}(n)\).
Appendix: Vectorization
With a bit of linear algebra, we can compute our recursive sum quite efficiently. Here is how:
Suppose we represent the coefficients of the expression \(S_p{n}\) as a matrix \(M\) such that \((i, j)^\text{th}\) entry of \(M\) is the coefficient of \(n^i\) in the expression for \(S_{j}(n)\). Then if the matrix is populated for all the coefficients of \(S_p(n)\) for \(p=1, \dots, k\), we can compute the \(k+1\) column using one matrix to vector multiplication and some vector additions.
The L.H.S. of the formula consists of two parts: \(nS_k(n)\) and \(-(a_{k,2} - 1)S_k(n) - \sum_{j=3}^{k+1} a_{k,j}S_{k+2-j}(n)\). The first part is just the formula for \(S_k(n)\) with all the exponents of the \(n\) terms incremented by 1. The second part can be computed using a matrix-vector multiplication.
I will illustrate this using a \(3\times 3\) matrix. Initially, it is occupied by
the coefficients of \(S_1(n)\) only.
\[
M = \begin{pmatrix}
\frac{1}{2} & 0 & 0\\\
\frac{1}{2} & 0 & 0\\\
0 & 0 & 0\\\
0 & 0 & 0
\end{pmatrix}
\]
For getting the second column, we will populate it first by \(nS_1(n)\). This is just the coefficients of the first column with one \(0\) added to the beginning.
So our matrix becomes
\begin{pmatrix}
\frac{1}{2} & 0 & 0\\\
\frac{1}{2} & \frac{1}{2} & 0\\\
0 & \frac{1}{2} & 0\\\
0 & 0 & 0
\end{pmatrix}
Now we want to find out the coefficients we need to subtract from the second column. For that we consider the column vector
\begin{pmatrix}
a_{k,k+1} \\\
\vdots\\\
a_{k,3}\\\
a_{k,2} - 1\\\
\end{pmatrix}
This is just the \(k^\text{th}\) column with the \(1\) subtracted from the last element.
For \(k=1\), this is just the one-element matrix \((-\frac{1}{2})\).
We multiply the upper-left \(2\times 1\) submatrix with the column matrix to get
\begin{pmatrix}
-\frac{1}{4}\\\
-\frac{1}{4}
\end{pmatrix}
This we need to subtract from upper \(2 \times 1\) sub-matrix of the second column to get
\begin{pmatrix}
\frac{1}{2} & \frac{1}{4} & 0\\\
\frac{1}{2} & \frac{3}{4} & 0\\\
0 & \frac{1}{2} & 0\\\
0 & 0 & 0
\end{pmatrix}
So the second column now corresponds to \((a_{1,1}+1)S_2(n)\). To get the column for \(S_2(n)\), we need to divide the whole column by \((a_{1,1}+1) = \frac{3}{2}\). The final matrix is
\begin{pmatrix}
\frac{1}{2} & \frac{1}{6} & 0\\\
\frac{1}{2} & \frac{1}{2} & 0\\\
0 & \frac{1}{3} & 0\\\
0 & 0 & 0
\end{pmatrix}
This means that \(S_2(n) = \frac{n}{6} + \frac{n^2}{2} + \frac{n^3}{3}\).
To compute the last column, we proceed in a similar manner. Duplicating the second column onto the third column and shifting it down by one row, we get
\begin{pmatrix}
\frac{1}{2} & \frac{1}{6} & 0\\\
\frac{1}{2} & \frac{1}{2} & \frac{1}{6}\\\
0 & \frac{1}{3} & \frac{1}{2}\\\
0 & 0 & \frac{1}{3}
\end{pmatrix}
Now, we need to multiply the \(3 \times 2\) upper left sub-matrix with the \(2 \times 1\) column matrix
\begin{align}
\begin{pmatrix}
\frac{1}{6}\\\
\frac{1}{2} - 1
\end{pmatrix} &= \begin{pmatrix}
\frac{1}{6}\\\
-\frac{1}{2}
\end{pmatrix}
\end{align}
Result is the \(3 \times 1\) column matrix
\begin{pmatrix}
0\\\
-\frac{1}{6}\\\
-\frac{1}{6}
\end{pmatrix}
Subtracting from the third column, we get \(M\) as
\begin{pmatrix}
\frac{1}{2} & \frac{1}{6} & 0\\\
\frac{1}{2} & \frac{1}{2} & \frac{1}{3}\\\
0 & \frac{1}{3} & \frac{2}{3}\\\
0 & 0 & \frac{1}{3}
\end{pmatrix}
Dividing the column by \((a_{2,1}+1) = \frac{4}{3}\), we get the final matrix
\begin{pmatrix}
\frac{1}{2} & \frac{1}{6} & 0\\\
\frac{1}{2} & \frac{1}{2} & \frac{1}{4}\\\
0 & \frac{1}{3} & \frac{1}{2}\\\
0 & 0 & \frac{1}{4}
\end{pmatrix}
This means that \(S_3(n) = \frac{n^2}{4} + \frac{n^3}{2} + \frac{n^4}{4}\).