We define the following polynomials, for $n≥0$: $$P_n(x)=(x+1)^{n+1}-x^{n+1}=\sum_{k=0}^{n}{\binom{n+1}{k}x^k}$$ For $n=0,1,2,3$ this gives us, $$P_0(x)=1\enspace P_1(x)=2x+1\enspace P_2(x)=3x^2+3x+1\enspace P_3(x)=4x^3+6x^2+4x+1$$
We then define the set $P_{(3)}=\{P_0,P_1,P_2,P_3\}$. It can be easily shown that this set is a basis over the vector space of polynomials of degree $3$ and lower. We take $3$ for the sake of brevity.
Taking the coefficients of these polynomials and turning them into column vectors, we can construct the matrix (coefficients from the lowest term to the highest term) $$\large{M_{P_{(3)}}}=\begin{pmatrix} 1 & 1 & 1 & 1 \\ 0 & 2 & 3 & 4 \\ 0 & 0 & 3 & 6 \\ 0 & 0 & 0 & 4 \end{pmatrix}$$ We'll call this matrix the pascal in the context of this post, and the above polynomials as pascal polynomials. The inverse of this matrix is the matrix, $$M_{P_{(3)}}^{-1}=\begin{pmatrix} 1 & -\frac{1}{2} & \frac{1}{6} & 0 \\ 0 & \frac{1}{2} & -\frac{1}{2} & \frac{1}{4} \\ 0 & 0 & \frac{1}{3} & -\frac{1}{2} \\ 0 & 0 & 0 & \frac{1}{4} \end{pmatrix}$$
We'll factor this matrix into two matrices as follows:
$$M_{P_{(3)}}^{-1}=\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & \frac{1}{2} & 0 & 0 \\ 0 & 0 & \frac{1}{3} & 0 \\ 0 & 0 & 0 & \frac{1}{4} \end{pmatrix}×\begin{pmatrix} 1 & -\frac{1}{2} & \frac{1}{6} & 0 \\ 0 & 1 & -1 & \frac{1}{2} \\ 0 & 0 & 1 & -\frac{3}{2} \\ 0 & 0 & 0 & 1 \end{pmatrix}$$ We can see the Bernoulli numbers in the first row of the matrix. Every column is a coefficient vector of a Bernoulli polynomial
The following are extended versions of these matrices:
Are there accepted names for these matrices and polynomials? What is the meaning of these relationships?
In particular, is there some treatment of using these matrices as change of basis transformations between representations of polynomials? E.g. from a linear combination of pascal polynomials to a linear combination of monomial terms.
Here's another way to look at this. The Bernoulli polynomials can be defined by the property
$$\int_x^{x+1} B_n(u) \, du = x^n.$$
So if we let $T$ be the operator from the set of polynomials to itself given by $(Tf)(x) = \int_x^{x+1} f(u) \, du$, then we have $(TB_n)(x) = x^n$. The operator $T$ sends $x^n$ to $$\int_x^{x+1} u^n \, du = \frac{1}{n+1}\left((x+1)^{n+1} - x^{n+1}\right) = \sum_{k=0}^{n} \frac{1}{n+1} {n+1 \choose k}x^k.$$
Writing $T$ as an infinite matrix with respect to the basis $1,x,x^2, \ldots$, gives
$$T = \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrr} 1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \ldots\\ 0 & 1 & 1 & 1 & \ldots\\ 0 & 0 & 1 & \frac{3}{2} & \ldots\\ 0 & 0 & 0 & 1& \ldots\\ \vdots & \vdots & \vdots & \vdots & \ddots \end{array}\right). $$
We can approximate the inverse by taking the inverse of this truncation, giving the matrix form of $T^{-1}$:
$$T^{-1} = \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrr} 1 & -\frac{1}{2} & \frac{1}{6} & 0 & \ldots\\ 0 & 1 & -1 & \frac{1}{2} & \ldots\\ 0 & 0 & 1 & -\frac{3}{2} & \ldots\\ 0 & 0 & 0 & 1& \ldots\\ \vdots & \vdots & \vdots & \vdots & \ddots \end{array}\right)$$
This operator sends $x^n$ to $B_n(x)$, so the columns are the coefficients of Bernoulli polynomials.
To see how this fits with your equations, note that we may factor the first matrix $T$ to remove the fraction $\frac{1}{n+1}$ in the formula for its coefficients, since multiplying by a diagonal matrix on the left scales the columns of a matrix by the diagonal entries. This gives $$ T = \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrr} 1 & 1 & 1 & 1 & \ldots\\ 0 & 2 & 3 & 4 & \ldots\\ 0 & 0 & 3 & 6 & \ldots\\ 0 & 0 & 0 & 4 & \ldots\\ \vdots & \vdots & \vdots & \vdots & \ddots \end{array}\right)\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrr} 1 & 0 & 0 & 0 & \ldots\\ 0 & 1/2 & 0 & 0 & \ldots\\ 0 & 0 & 1/3 & 0 & \ldots\\ 0 & 0 & 0 & 1/4 & \ldots\\ \vdots & \vdots & \vdots & \vdots & \ddots \end{array}\right)$$
The matrix on the left is your $M_P$ (infinitely extended). Calling the diagonal matrix on the right $D$, we have $T = M_pD$, or $DT^{-1} = M_P^{-1}$ which is your last equation.
To my mind, this is the most natural way to compute Bernoulli polynomials. A few years ago I was playing around with this idea before I knew what a Bernoulli polynomial was. I was slightly disappointed, although not too surprised, to hear that someone else had discovered this first. I was beaten by some three hundred years, no less. It ended up as part of my undergraduate thesis. :)