A relationship between matrices, bernoulli polynomials, and binomial coefficients 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:  


*

*Pascal Matrix

*Inverse Pascal Matrix

*Last, Bernoulli Matrix
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.
 A: 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. :)
A: This is a very fascinating problem, but in the $P_{(4)}$ case, we have
$$\tag{1}
M_{P_{(4)}}^{-1}=\begin{pmatrix}
1 & -1/2 &  1/6 &   0  & -1/30\\
  &  1/2 & -1/2 &  1/4 &   0\\
  &      &  1/3 & -1/2 &  1/3\\
  &      &      &  1/4 & -1/2\\
  &      &      &      &  1/5 
\end{pmatrix}
$$
This factors out as
$$
M_{P_{(4)}}^{-1}=\begin{pmatrix}
1 & -1 &  1 &   0  & -1/6\\
  &  1 & -3/2 &  1 &   0\\
  &      &  1 & -2 &  5/3\\
  &      &      &  1 & -5/2\\
  &      &      &      &  1
\end{pmatrix}
\begin{pmatrix}
1 &  &   &     & \\
  &  1/2 &  &   &   \\
  &      &  1/3 & &  \\
  &      &      &  1/4 & \\
  &      &      &      &  1/5 
\end{pmatrix}
$$
Unfortunately, as we can read off, the columns in the matrix on the left do not give us Bernoulli polynomials :(
Although, please note the first row of the matrix in (1) gives us Bernoulli numbers. Perhaps the first row of the inverse matrix does give you the Bernoulli numbers, but you do not obtain the Bernoulli polynomials.
Again, read this with caution and suspicion: the good reader should compute this and verify it for him or herself!
Addendum
Write out
$$
M_{P_{(n)}} = DU
$$
where $D$ is a diagonal matrix, and $U=I+X$ is a unit upper triangular matrix (here $X$ is the strictly upper triangular part of $U$). The claim is that
$$(DU)^{-1}=(I+X)^{-1}D^{-1}=(I-X+X^{2}-X^{3}+\dots+(-1)^{n}X^{n})D^{-1}$$
produces the coefficients of Bernoulli polynomials in the columns, and the first row is the Bernoulli numbers. I claim the first (Bernoulli polynomial data) implies the latter (the first row consists of Bernoulli numbers), and claim this is obvious.
Conjecture
It seems that the OP found the matrix version for the Bernoulli polynomials described on Wikipedia using forward differences.
I would have urged the OP to change notation, and use lower triangular matrices. Why? Because then you could write
$$
M_{P_{(n)}}\begin{bmatrix}0\\ 0\\\vdots\\0\\ x^{n}\end{bmatrix} \cong (1+x)^{n+1}-x^{n+1}
$$
in vector form. Current notation demands we use row vectors for polynomials.
I'm about to go to bed, so I do not believe I have time to prove my conjecture. If no one has proven it by tomorrow, I'll try writing up a proof.
A: (Just stumbled on this question.)
These matrix relationships are easily understood in terms of the lower triangular coefficient matrices and the umbral calculus of Sheffer polynomials.
The Bernoulli polynomials are an Appell sequence of Sheffer polynomials with the e.g.f., in umbral notation,
$$e^{B.(x)t} = B(t) e^{xt} = \frac{t}{e^t-1} \; e^{xt} = e^{b.t} \; e^{xt} = e^{(b.+x)t}$$
with $(b.)^n = b_n$, the Bernoulli numbers, giving the Bernoulli polynomials
$$(B.(x))^n = B_n(x) = (b.+x)^n = \sum_{k=0}^n \; \binom{n}{k} \; b_{n-k} \; x^k$$
$$ = B(D)\; x^n = e^{b.D} \; x^n = \frac{D}{e^D-1} \; x^n,$$
where $D = d/dx$,
and the associated lower triangular matrix $[B]$ with elements
$$B_{n,k} = \binom{n}{k} \; b_{n-k}.$$
Define the reciprocal polynomials $R_n(x)$ with the e.g.f.
$$e^{R.(x)t} = R(t) e^{xt} = \frac{1}{B(t)} \; e^{xt}  =  \frac{e^t-1}{t} \; e^{xt} = e^{r.t} \; e^{xt} = e^{(r.+x)t}$$
with $r_n = \frac{1}{n+1}$, giving the reciprocal polynomials
$$(R.(x))^n = R_n(x) = (r.+x)^n = \sum_{k=0}^n \; \binom{n}{k} \; \frac{1}{n-k+1} \; x^k$$
$$ = R(D)\; x^n = e^{r.D} \; x^n = \frac{e^D-1}{D} \; x^n$$
$$=\frac{(1+x)^{n+1}-x^{n+1}}{n+1}= \frac{P_n(x)}{n+1}$$
and the associated lower triangular matrix $[R]$ with elements
$$R_{n,k} = \binom{n}{k} \; r_{n-k} = \binom{n}{k} \; \frac{1}{n-k+1} .$$
The pair of Bernoulli and reciprocal polynomial sequences are an umbral inverse pair under umbral composition; that is,
$$B_n(R.(x)) = x^n = R_n(B.(x))$$
$$= \sum_{k=0}^n B_{n,k} \; R_k(x) = \sum_{k=0}^n B_{n,k} \; \sum_{j=0}^k R_{k,j} \; x^j$$
$$= \sum_{k=0}^n R_{n,k} \; B_k(x) = \sum_{k=0}^n R_{n,k} \; \sum_{j=0}^k B_{k,j} \; x^j.$$
These relations imply the pair of matrices $[R]$ and $[B]$ are an inverse pair; i.e.,
$$[R]^{-1} = [B],$$
so inverting $[R]$ gives the Bernoulli numbers in the first column of $[R]^{-1}$.
To simply show the umbral inverse relation holds, note
$$B_n(R.(x)) = R(D) B_n(x)= R(D) B(D) x^n = R(D) \frac{1}{R(D)}x^n = x^n$$
and similarly
$$ R_n(B.(x)) = B(D) R(D) x^n = x^n.$$
These machinations easily generalize to other pairs of Sheffer polynomials and their associated lower triangular coefficient matrices, such as the Stirling polynomials of the first and second kinds, so important in combinatorics and analysis.
See OEIS A135278, A074909, and A133314 for related matrices, analysis, and refs.
A: Regarding the interpretation of the matrices: 
Let 
$$X = 
\left(\begin{array}{c}
1 \\ 
x \\
\vdots \\
x^n 
\end{array}\right) 
\qquad \textrm{and}\qquad  
A =
\left(\begin{array}{c}
a_0 \\ 
a_1 \\ 
\vdots \\ 
a_n
\end{array}\right)$$ 
so 
$$X^T A = \sum_{k=0}^n a_k x^k.$$
(I've written $X^T A$ instead of the preferred $A^T X$ to agree with your notation and avoid things like $(M^{-1})^T$.)
Then 
$$\begin{eqnarray*}
X^T A &=& X^T M M^{-1} A \\
&=& {X'}^T A'
\end{eqnarray*}$$
where 
$X' = M^T X$
is the new basis and $A' = M^{-1} A$ are the coefficients in the new basis.
This interpretation holds for any invertible transformation $M$.
The formalism can be generalized in a straightforward way to transformations between any two bases of polynomials.
For your transformation we find the coefficients of 
$\sum_{k=0}^3 a_k x^k$ in the "Pascal basis" are 
$$M_{P_{(3)}}^{-1}A = \left(
\begin{array}{c}
a_0 - \frac{1}{2} a_1 + \frac{1}{6} a_2 \\
\frac{1}{2}a_1 - \frac{1}{2} a_2 + \frac{1}{4} a_3 \\
\frac{1}{3} a_2 - \frac{1}{2} a_3 \\
\frac{1}{4} a_3
\end{array}
\right).$$
Thus, 
$$\begin{eqnarray*}
\sum_{k=0}^3 a_k x^k &=& 
\left(a_0 - \frac{1}{2} a_1 + \frac{1}{6} a_2\right) P_0(x) \\
&&+ \left(\frac{1}{3} a_2 - \frac{1}{2} a_3\right) P_1(x) \\
&&+ \left(\frac{1}{3} a_2 - \frac{1}{2} a_3\right) P_2(x) \\
&&+ \frac{1}{4} a_3 P_3(x).
\end{eqnarray*}$$
A: Given a polynomial $q(x)$ of degree $n+1$, the derivative, $Dq(x)=\frac{d}{dx}q(x)$, and the forward difference, $\Delta q(x)=q(x+1)-q(x)$, are both polynomials of degree $n$. Elementary calculus describes the easy manipulation needed to recover $q(x)$ from $Dq(x)$, up to a constant. By linear algebra—essentially the calculations you have done—it is possible to recover $q(x)$ from $\Delta q(x)$, up to a constant.
These parallels strongly suggest that there is a systematic relationship between $\Delta q(x)$ and $Dq(x)$, and that one should focus on the transformation that turns one into the other. If $\Delta q(x)$ can be turned into $Dq(x)$, then $q(x)$ can be recovered immediately, up to a constant, by taking the antiderivative. The needed transformation is, in fact, given by your Bernoulli matrix, extended arbitrarily, which we may call $B$, and provides as good a definition of the Bernoulli polynomials as any. If $G$ is the vector of coefficients of $\Delta q(x)$, in the monomial basis, and $F'$ the vector of coefficients of $Dq(x)$, in the same basis, then
$$
F'=BG.
$$
The effect of the linear operator $B$ is to replace each monomial $x^k$ with the corresponding Bernoulli polynomial $B_k(x)$. The operators $D$ and $\Delta$ may be represented in this basis by the semi-infinite matrices
$$
D=\begin{bmatrix}0 & 1 & 0 & 0 & 0 & \cdots\\ 0 & 0 & 2 & 0 & 0 & \cdots\\ 0 & 0 & 0 & 3 & 0 & \cdots\\ 0 & 0 & 0 & 0 & 4 & \cdots\\ \vdots & \vdots & \vdots & \vdots & \ddots & \ddots\end{bmatrix}\quad\text{and}\quad\Delta=\begin{bmatrix}0 & 1 & 1 & 1 & 1 & \cdots\\ 0 & 0 & 2 & 3 & 4 & \cdots\\ 0 & 0 & 0 & 3 & 6 & \cdots\\ 0 & 0 & 0 & 0 & 4 & \cdots\\ \vdots & \vdots & \vdots & \vdots & \ddots & \ddots\end{bmatrix}.
$$
With $F$ as the vector of coefficients of $q(x)$ in the monomial basis, we have
$$
DF=B\Delta F.
$$
The operators $D$ and $\Delta$ are singular, while $B$ is invertible. The null spaces of $D$ and $\Delta$ are the same, as they must be, and correspond to the set of constant polynomials. Since $F$ is arbitrary, we have $D=B\Delta$. As a consequence of linearity and the chain rule, the operators $D$ and $\Delta$ commute.
Applying $D=B\Delta$ to $F'$, we have $D^2F=B\Delta DF=BD\Delta F$. But $D^2F=DB\Delta F$, so we have $DBG=BDG$ for arbitrary $G$—every $G$ is the forward difference of some $F$—and we conclude that $D$ and $B$ commute. This is essentially the property
$$
B'_n(x)=nB_{n-1}(x).
$$
Using $D=B\Delta$ in $DB=BD$, we have $B\Delta B=B^2\Delta$, from which it follows, since $B$ is invertible, that $B$ and $\Delta$ commute. This is essentially the property that
$$
\Delta Bx^n=B_n(x+1)-B_n(x)=B\Delta x^n=Dx^n=nx^{n-1},
$$
and also implies
$$
DF=\Delta BF.
$$
I'll stop the development here, but the three mutually commuting operators $D$, $\Delta$, and $B$ can be used to find the generating function for the Bernoulli polynomials and to derive a great many other properties.
Some details: To address the questions in your final two paragraphs: I am not aware of any standard names for the three matrices defined in the post; the term "Bernoulli matrix" is sometimes used for a completely unrelated concept. You will recognize $M_{P_{(n)}}$ as a submatrix of $\Delta$, and the factor on the left in your factorization of $M_{P_{(n)}}^{-1}$ as the inverse of a submatrix of $D$.

*

*The Pascal matrix $M_{P_{(n)}}$ acts on the space of polynomials of degree $n+1$ with zero constant term, or, better, on the space of equivalence classes of polynomials of degree $n+1$, where two polynomials are considered equivalent when they differ by a constant. It produces the polynomial of degree $n$ equal to the forward difference of any equivalence class representative. Let $f(x)=f_{n+1}x^{n+1}+f_nx^n+\ldots+f_0$ be any representative of an equivalence class and define the column vector $F$ by $F^T=\begin{bmatrix}f_1 & f_2 & \ldots & f_{n+1}\end{bmatrix}$. Let $g(x)=\Delta f(x)=f(x+1)-f(x)=g_nx^n+g_{n-1}x^{n-1}+\ldots+g_0$ and define the column vector $G$ by $G^T=\begin{bmatrix}g_0 & g_1 & \ldots & g_n\end{bmatrix}$. Then
$$
G=M_{P_{(n)}}F.
$$

*The inverse Pascal matrix reverses this process, producing the equivalence class of polynomials that have a given forward difference:
$$
F=M_{P_{(n)}}^{-1}G.
$$

*You have factorized the inverse Pascal matrix as $M_{P_{(n)}}^{-1}=AB$, where $A=\operatorname{diag}(\begin{bmatrix}1 & \frac{1}{2} & \ldots & \frac{1}{n+1}\end{bmatrix})$ and $B$ is your Bernoulli matrix. Considering the Bernoulli matrix first, it acts on the space of polynomials of degree $n$ and produces a polynomial in the same space. What it does is to transform the forward difference of a polynomial into the derivative of that polynomial:
$$
F'=BG,
$$
where ${F'}^T=\begin{bmatrix}f_1 & 2f_2 & \ldots & (n+1)f_{n+1}\end{bmatrix}$ is the vector of coefficients of $f'(x)$. It accomplishes this by replacing $x^k$ with $B_k(x)$. The matrix $A$ takes a polynomial of degree $n$ and produces its antiderivative, which lives in the set of equivalence classes of polynomials of degree $n+1$.

Point 1 above follows from your definition of $M_{P_{(n)}}$ and linearity of the forward difference operator. Point 2 is obvious; invertibility is clear because $M_{P_{(n)}}$ is upper-triangular with nonzero elements on the diagonal. To compute $M_{P_{(n)}}^{-1}$ and establish the factorization in Point 3, let the first row (row $0$) of $M_{P_{(n)}}^{-1}$ be $\begin{bmatrix}b_0 & b_1 & b_2 & \ldots\end{bmatrix}$. Forming the product of this row with column $j$ of $M_{P_{(n)}}$ gives
$$
\sum_{i=0}^j\binom{j+1}{i}b_i=\delta_{j,0},
$$
a well-known recurrence for the Bernoulli numbers. More generally, the first $k$ elements of row $k$ of $M_{P_{(n)}}^{-1}$ are zero, so we may write this row as $\begin{bmatrix}0 & \ldots & 0 & b_0^{(k)} & b_1^{(k)} & \ldots\end{bmatrix}$. The product of this row with column $j$ of $M_{P_{(n)}}$ is $0$ for $j<k$; for column $j+k$, $j\ge0$ the product is
$$
\delta_{k,j+k}=\sum_{i=0}^j\binom{j+k+1}{i+k}b_i^{(k)}=\binom{j+k+1}{k}\sum_{i=0}^j\binom{j+1}{i}\frac{b_i^{(k)}}{\binom{i+k}{k}},
$$
where the second equality follows from the identity
$$
\binom{a+c}{b+c}=\binom{a}{b}\frac{\binom{a+c}{c}}{\binom{b+c}{c}}.
$$
Therefore
$$
\sum_{i=0}^j\binom{j+1}{i}\frac{(k+1)b_i^{(k)}}{\binom{i+k}{k}}=\frac{(k+1)\delta_{k,j+k}}{\binom{j+k+1}{k}}=\delta_{j,0},
$$
and the quantities $(k+1)b_i^{(k)}/\binom{i+k}{k}$, $i\ge0$, satisfy the same system of equations as the quantities $b_i$ do. We conclude that
$$
b_i^{(k)}=\frac{1}{k+1}\binom{i+k}{k}b_i.
$$
The first factor in this expression accounts for the matrix $A$. Comparing the remaining factors with the standard expression for the Bernoulli polynomials in terms of binomial coefficients and Bernoulli numbers establishes that the columns of $B$ are the sequences of coefficients of Bernoulli polynomials.
Embedding $M_{P_{(n)}}$ and $A^{-1}$ in the matrices $\Delta$ and $D$, as in the introduction to this post, allows us to dispense with worrying about degrees of polynomials and equivalence classes, and frees us to work simply with the space of all polynomials.
