Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

First Variant:

Bernoulli numbers can easily be expressed by linear algebra equations. For example just using the recursion formula

$$\sum_{k=0}^{n-1}{n\choose k}B_k=0$$

which is equation (34) from the MathWorld page. Notice that the sum only goes to $n-1$.

We can formulate the calculation of the Bernoulli numbers as a null space calculation as follows:


$$\mathbf{L} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0\\ 1 & 1 & 0 & 0 & 0 & 0\\ 1 & 2 & 1 & 0 & 0 & 0\\ 1 & 3 & 3 & 1 & 0 & 0\\ 1 & 4 & 6 & 4 & 1 & 0\\ 1 & 5 &10 & 10& 5 & 1\end{bmatrix}$$

be the $6x6$ lower triangular Pascal matrix. Set $l_{21} = 0$, then you will find the vector $b=[1, -1/2, 1/6, 0, -1/30, 0]^T$ as a solution to $(L-I)b=0$. This vector is a vector of Bernoulli numbers.

Second Variant:

However there is another way to find the Bernoulli numbers with matrices. Consider

$$\mathbf{A} = \begin{bmatrix} -\frac{1}{2} & -\frac{1}{6} & -\frac{1}{12} & -\frac{1}{20} & -\frac{1}{30} & -\frac{1}{42}\\ \frac{1}{1} & 0 & 0 & 0 & 0 & 0\\ 0 & \frac{1}{2} & 0 & 0 & 0 & 0\\ 0 & 0 & \frac{1}{3} & 0 & 0 & 0\\ 0 & 0 & 0 & \frac{1}{4} & 0 & 0\\ 0 & 0 & 0 & 0 & \frac{1}{5} & 0 \end{bmatrix}$$

So this matrix was constructed by putting $-1/(1\cdot2), -1/(2\cdot3), -1/(3\cdot4), ...$ on the first row and $1/1, 1/2, 1/3, ...$ on the subdiagonal.

Now with $k!A^ke$ with $e:=[1, 0, 0, 0, 0, 0]^T$ and $k < 6$ we will get the coefficients of the kth Bernoulli polynomial. With the constant term in the 'first' (leftmost) position. For example

$$3!A^3e = [0, 1/2, -3/2, 1, 0, 0]^T$$


$$4!A^4e = [-1/30, 0, 1, -2, 1, 0]^T$$

Of course there is a difference between the coefficients of the Bernoulli polynomials and the Bernoulli numbers - but it isn't very big since

$$B_n(x) = \sum_{k=0}^n{n\choose k}B_kx^{n-k}$$

this is from wikipedia.

Clearly the inspiration for this slightly different way of obtaining the Bernoulli numbers/polynomials was the 'calculus' definition of the Bernoulli polynomials:

$$B'_k(w) = k \cdot B_{k-1}(w)$$

$$\int_0^1 B_k(t)\textrm{d}t = B_{k+1}(1)-B_{k+1}(0)=0$$

But how are these two methods related? Is there maybe a linear algebra way of showing that these calculations result in the same (bernoulli) vectors?

EDIT 1: The second variant (defined in this question) is consistent with $B(s)=-s\zeta (1-s)$ since you can just calculate $\Gamma (s+1)A^s e$ - and it gets closer for bigger $n$ (this is for real $s$ - complex $s$ is untested). I don't know how you could make the first variant consistent with $B(s)=-s\zeta (1-s)$ yet...

EDIT 2: Let $$D:= \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 2 & 0 & 0 & 0\\ 0 & 0 & 0 & 3 & 0 & 0\\ 0 & 0 & 0 & 0 & 4 & 0\\ 0 & 0 & 0 & 0 & 0 & 5\\ \end{bmatrix} $$

be the $6x6$ 'differentiation operator' for polynomials in monomial basis, where the polynomial vectors are given with the constant term at the top (and the highest degree power at the bottom). Then if $b_k=k!A^ke$ is a vector of coefficients of the kth bernoulli polynomial produced via the 'second variant' we have for example $Db_k=kb_{k-1}$.

But another relationship that the bernoulli polynomials obey is more useful here. We know that the forward difference of the nth bernoulli polynomial results in the derivative of the the monomial of power $n-1$ $$B_n(x+1)-B_n(x)=nx^{n-1}$$ from wikipedia. So if we change the basis of a polynomial from monomial to 'bernoulli' the forward difference operator becomes the derivative operator. But the forward difference operator for polynomials in monomial basis is just $e^D-I$ - where $e^D$ is the upper triangular pascal matrix again from wikipedia. The 'change of basis matrix' from monomial to 'bernoulli' is just $B:=[b_0,b_1,...,b_5]$ with $b_k=k!A^ke.$ So $$(e^D-I)B=D.$$ (btw matrix B has full rank $6$ - $(e^D-I)$ and $D$ only have rank $5$).

Since $B$ has full rank we can write $e^D-I = DB^{-1}$ and because $D$ is nilpotent the left side is a finite (matrix) sum, namely $$\frac{1}{1!}D^1+\frac{1}{2!}D^2+\frac{1}{3!}D^3+\frac{1}{4!}D^4+\frac{1}{5!}D^5=DB^{-1}$$ So if we multiply both sides by the 'integration operator' for polynomials we get $$B^{-1}=\sum_{k=0}^4\frac{D^k}{(k+1)!}.$$ However since the 'integration operator' for polynomials with all zeros on the top row doesn't have full rank, it seems that one value of $B^{-1}$ remains undetermined (although $B^{-1}$ calculated this way is still invertible and you can still get the first four bernoulli polynomials by inverting $B^{-1}$).

From this sum we can conlude that $B^{-1}-I$ is nilpotent so $(B^{-1}-I)^5=0$ or $$\sum_{k=0}^5 (-1)^{k}{5\choose k}(B^{-1})^k=0.$$

The first term in this sum is the identity, so we can solve for that and multiply both sides with $B$ to get

$$B = -\sum_{k=1}^5 (-1)^{k}{5\choose k}(B^{-1})^{k-1}.$$

This is getting very close to the exlicit formula we have for the bernoulli polynomials see wikipedia. However that formula works for the values of the bernoulli polynomials. The formula presented here just uses the coefficients. (Also it not only considers the coefficients, but the coefficients of the fist $n=5$ bernoulli polynomials all at once - which is awkward).

The formula $(e^D-I)B=D$ provides a slightly stronger connection between the two variants since the left side looks a lot like our modified lower triangular pascal matrix and the right side looks a lot like $A$ since $A$ is just the 'integral operator' for polynomials where the first row fixes the constant term (definite integration). Even though this doesn't completely answer the question, i think it is still worth noting.

PS: you should have a look at the article that Gottfried Helms linked in the comments - in it he introduces the diagonal matrix $J$ whose diagonal entries just alternate between $+1$ and $-1$ (starting with $+1$). Then he makes the observation that $JLb=b$ (where $L$ is the unaltered lower triangular Pascal matrix). So $b$ is in the null space of $JL-I$. Note that $JL$ is almost (but not quite) the inverse of $L$. This seems entirely compatible with the 'first variant' defined here. However the null space of $JL-I$ is a dimension bigger than the one of $L-I$ (with altered $L$) so the Bernoulli vector doesn't turn out to be a (possibly scaled) basis vector of the null space, but just in the null space somewhere.

share|cite|improve this question
As $L$ is a matrix with numerical entries, and not a function of two variables, I don't know what "Set $L(2,1)=0$" means. Nor do I understand how anything can be a solution to "$L-I=0$". – Gerry Myerson Jul 18 '11 at 0:37
@Gerry: My naïve interpretation of Peter is that the eigenvector of the Pascal matrix with the entry in the second row and first column set to 0 corresponding to the eigenvalue 1 can be normalized such that its components are Bernoulli numbers... – Jerry Jul 18 '11 at 5:49
Peter: the more conventional notation is that if you have some matrix $\mathbf L$, then you want $l_{21}=0$ (i.e. use doubly subscripted small letters for matrix entries if you're representing matrices as boldface capital letters...) – Jerry Jul 18 '11 at 5:51
@Peter, so, when you write $L-I=0$, you really mean $(L-I)b=0$, right? It would be a lot easier to come to grips with your question if you would write what you really mean, instead of sorta kinda something like what you mean. – Gerry Myerson Jul 18 '11 at 6:40
@PeterSheldrick : A slightly extended and generalized discussion of the matrix-representation of Bernoulli-polynomials and the diagonalization-aspect is in… which I forgot when I answered in July - maybe you'll find some more in it. – Gottfried Helms Oct 4 '11 at 8:32
up vote 1 down vote accepted

Perhaps the thoughts which I followed recently give a satisfactory answer. I've looked at the "ZETA"-matrix and fiddled out the connection to the integral-representation in the Euler-MacLaurin-Formula. The article is not yet ready, needs some brushing and completing references and such. But as it might be helpful here: here is the link the integral in Euler-MacLaurin in connection with the Pascal-matrix

[update]: added the re-translation into bernoulli-numbers in the final formulae in the text

share|cite|improve this answer
Just "normalizing" - I'm studying zeta-references, indefinite summation, iteration and series in various contexts; to find identities it is useful to reduce things to "normalized" expressions using the same set of (possibly fewest) symbols. Here I prefer $\zeta$ over bernoulli-numbers because they also allow generalizations to fractional indexes where this is meaningful. But for the casual reader I'll surely add the expressions in terms of Bernoulli-numbers in the article, too. – Gottfried Helms Jul 20 '11 at 14:11

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.