A Simple Method to Compute the Sum of Powers of Natural Numbers

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 )

Assuming that we know that the formula is a polynomial of degree (which happens to be the case), we can easily derive the coefficients by substituting for , 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 as .

Case

We will take our base case as . Obviously, is just .

Case

To compute the formula for , consider the following sum:

Now, is the same as . But as there are number of , number of , etc., we can write also in the following manner:

Now we can equate the two expressions for , and get:

Which checks out!

Case

Rewriting as in the previous case, we see that

Equating, we can obtain the formula for .

General Case

Now we can generalize what we have been doing so far. But first, we need to assume that is a polynomial of degree 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 is a polynomial of degree , we will prove that is similarly a polynomial of degree 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 , for all . We can take one of as the base case.

For computing the formula for , consider the expression:

But also,

Equating both the expressions, we get

From the last expression we can see that the largest power of is times the largest power of . Since is a polynomial of degree , must be a polynomial of degree . Also note that all the coefficients of the recursive sum are just the coefficients of . So if we have , we can readily compute the expression for .

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 as a matrix such that entry of is the coefficient of in the expression for . Then if the matrix is populated for all the coefficients of for , we can compute the column using one matrix to vector multiplication and some vector additions.

The L.H.S. of the formula consists of two parts: and . The first part is just the formula for with all the exponents of the terms incremented by 1. The second part can be computed using a matrix-vector multiplication.

I will illustrate this using a matrix. Initially, it is occupied by the coefficients of only.

For getting the second column, we will populate it first by . This is just the coefficients of the first column with one added to the beginning.

So our matrix becomes

Now we want to find out the coefficients we need to subtract from the second column. For that we consider the column vector

This is just the column with the subtracted from the last element.

For , this is just the one-element matrix .

We multiply the upper-left submatrix with the column matrix to get

This we need to subtract from upper sub-matrix of the second column to get

So the second column now corresponds to . To get the column for , we need to divide the whole column by . The final matrix is

This means that .

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

Now, we need to multiply the upper left sub-matrix with the column matrix

Result is the column matrix

Subtracting from the third column, we get as

Dividing the column by , we get the final matrix

This means that .