We have,
$$\sum_{k=1}^n k^3 = \Big(\sum_{k=1}^n k\Big)^2$$
$$2\sum_{k=1}^n k^5 = -\Big(\sum_{k=1}^n k\Big)^2+3\Big(\sum_{k=1}^n k^2\Big)^2$$
$$2\sum_{k=1}^n k^7 = \Big(\sum_{k=1}^n k\Big)^2-3\Big(\sum_{k=1}^n k^2\Big)^2+4\Big(\sum_{k=1}^n k^3\Big)^2$$
and so on (apparently).
Is it true that the sum of consecutive odd $m$ powers, for $m>1$, can be expressed as sums of squares of sums* in a manner similar to the above? What is the general formula?
*(Edited re Lord Soth's and anon's comment.)
This is a partial answer, it just establishes the existence.
We have $$s_m(n) = \sum_{k=1}^n k^m = \frac{1}{m+1} \left(\operatorname{B}_{m+1}(n+1)-\operatorname{B}_{m+1}(1)\right)$$ where $\operatorname{B}_m(x)$ denotes the monic Bernoulli polynomial of degree $m$, which has the following useful properties: $$\begin{align} \int_x^{x+1}\operatorname{B}_m(t)\,\mathrm{d}t &= x^m \quad\text{(from which everything else follows)} \\\operatorname{B}_{m+1}'(x) &= (m+1)\operatorname{B}_m(x) \\\operatorname{B}_m\left(x+\frac{1}{2}\right) &\begin{cases} \text{is even in $x$} & \text{for even $m$} \\ \text{is odd in $x$} & \text{for odd $m$} \end{cases} \\\operatorname{B}_m(0) = \operatorname{B}_m(1) &= 0 \quad\text{for odd $m\geq3$} \end{align}$$
Therefore, $$\begin{align} s_m(n) &\text{has degree $m+1$ in $n$} \\s_m(0) &= 0 \\s_m'(0) &= \operatorname{B}_m(1) = 0\quad\text{for odd $m\geq3$} \\&\quad\text{(This makes $n=0$ a double zero of $s_m(n)$ for odd $m\geq3$)} \\s_m\left(x-\frac{1}{2}\right) &\begin{cases} \text{is even in $x$} & \text{for odd $m$} \\ \text{is odd in $x$} & \text{for even $m\geq2$} \end{cases} \end{align}$$
Consider the vector space $V_m$ of univariate polynomials $\in\mathbb{Q}[x]$ with degree not exceeding $2m+2$, that are even in $x$ and have a double zero at $x=\frac{1}{2}$. Thus $V_m$ has dimension $m$ and is clearly spanned by $$\left\{s_j^2\left(x-\frac{1}{2}\right)\mid j=1,\ldots,m\right\}$$ For $m>0$, we find that $s_{2m+1}(x-\frac{1}{2})$ has all the properties required for membership in $V_m$. Substituting $x-\frac{1}{2}=n$, we conclude that there exists a representation $$s_{2m+1}(n) = \sum_{j=1}^m a_{m,j}\,s_j^2(n) \quad\text{for $m>0$ with $a_{m,j}\in\mathbb{Q}$}$$ of $s_{2m+1}(n)$ as a linear combination of squares of sums.
Not an answer, but only some results I found.
Denote $s_m=\sum_{k=1}^nk^m.$ Suppose that $$s_{2m+1}=\sum_{i=1}^m a_is_i^2,a_i\in\mathbb Q.$$
Use the method of undetermined coefficients, we can get a list of $\{m,\{a_i\}\}$: $$\begin{array}{l} \{1,\{1\}\} \\ \left\{2,\left\{-\frac{1}{2},\frac{3}{2}\right\}\right\} \\ \left\{3,\left\{\frac{1}{2},-\frac{3}{2},2\right\}\right\} \\ \left\{4,\left\{-\frac{11}{12},\frac{11}{4},-\frac{10}{3},\frac{5}{2}\right\}\right\} \\ \left\{5,\left\{\frac{61}{24},-\frac{61}{8},\frac{28}{3},-\frac{25}{4},3\right\}\right\} \\ \left\{6,\left\{-\frac{1447}{144},\frac{1447}{48},-\frac{332}{9},\frac{595}{24},-\frac{21}{2},\frac{7}{2}\right\}\right\} \\ \left\{7,\left\{\frac{5771}{108},-\frac{5771}{36},\frac{5296}{27},-\frac{2375}{18},56,-\frac{49}{3},4\right\}\right\} \\ \left\{8,\left\{-\frac{53017}{144},\frac{53017}{48},-\frac{60817}{45},\frac{21817}{24},-\frac{3861}{10},\frac{1127}{10},-24,\frac{9}{2}\right\}\right\} \\ \left\{9,\left\{\frac{2755645}{864},-\frac{2755645}{288},\frac{632213}{54},-\frac{1133965}{144},\frac{13379}{4},-\frac{11725}{12},208,-\frac{135}{4},5\right\}\right\}\end{array}$$
Through some observation, we can find that $a_i a_{i+1}<0,a_m=\dfrac{m+1}2,a_2=-3a_1,a_{m-1}=\dfrac{1}{24} m^2 (m+1).$
Adding to @Hecke's answer
First step was to build the matrix $A$ of coefficients of the squares of the Faulhaber-polynomials. Next step, to build the matrix $B$ of the coefficients of the un-squared Faulhaber-polynomials, but only of the odd orders.
Then the matrix $C$ of coefficients for the composition, such that $C \cdot A = B$ occurs by inversion of $A$ such that $C = B \cdot A^{-1}$. This can be done with any truncation of size. Important: Note that the first two columns of $A$ and $B$ are completely zero, so we must reduce that matrices to the remaining columns.
The inverse of the coefficients-matrix $C$ (which itself can be seen in the inner braces of Hecke's answer) might also be interesting; at least it has a slightly simpler structure:
$$ C^{-1}= \tiny \begin{bmatrix}
2 & . & . & . & . & . & . & . & . & . & . & . & . & . \\
. & 1 & . & . & . & . & . & . & . & . & . & . & . & . \\
. & 1/3 & 2/3 & . & . & . & . & . & . & . & . & . & . & . \\
. & . & 1/2 & 1/2 & . & . & . & . & . & . & . & . & . & . \\
. & . & -1/15 & 2/3 & 2/5 & . & . & . & . & . & . & . & . & . \\
. & . & . & -1/6 & 5/6 & 1/3 & . & . & . & . & . & . & . & . \\
. & . & . & 1/21 & -1/3 & 1 & 2/7 & . & . & . & . & . & . & . \\
. & . & . & . & 1/6 & -7/12 & 7/6 & 1/4 & . & . & . & . & . & . \\
. & . & . & . & -1/15 & 4/9 & -14/15 & 4/3 & 2/9 & . & . & . & . & . \\
. & . & . & . & . & -3/10 & 1 & -7/5 & 3/2 & 1/5 & . & . & . & . \\
. & . & . & . & . & 5/33 & -1 & 2 & -2 & 5/3 & 2/11 & . & . & . \\
. & . & . & . & . & . & 5/6 & -11/4 & 11/3 & -11/4 & 11/6 & 1/6 & . & . \\
. & . & . & . & . & . & -691/1365 & 10/3 & -33/5 & 44/7 & -11/3 & 2 & 2/13 & . \\
. & . & . & . & . & . & . & -691/210 & 65/6 & -143/10 & 143/14 & -143/30 & 13/6 & 1/7
\end{bmatrix}$$
... and the numerator of $691$ in rows $13,14,15,\ldots$ surely rings some bell concerning Bernoulli-numbers
Note: this matrix gives the compositions of the squared sums-of-powers by the higher odd orders of the original sums-of-powers, so helps for the inverse of your problem:
$$ \begin{array}{rrrrr} S(0,x)^2 &=& 2 S(1,x) & - 1 S(0,x) &\\
\hline \\
S(1,x)^2 &=& 1 S(3,x) \\
S(2,x)^2 &=& 1/3 \; S(3,x) &+2/3 \; S(5,x) \\
S(3,x)^2 &=&& 1/2 \; S(5,x) &+1/2 \; S(7,x) \\
S(4,x)^2 &=&& -1/15 \; S(5,x)&+2/3 \; S(7,x)&+2/5 \; S(9,x) \\
\vdots \end{array} $$
The first row does not fit into the pattern. The coefficients are compositions of binomials and Bernoulli-numbers, which might best be written as below according to the analysis of the matrix $E$
[update] The following matrix displays coefficients after the Bernoulli-numbers and the column-index are removed. $$ D = \Tiny \begin{bmatrix} 2 & . & . & . & . & . & . & . & . & . & . & . \\ . & 2 & . & . & . & . & . & . & . & . & . & . \\ . & 4 & 2 & . & . & . & . & . & . & . & . & . \\ . & . & 9 & 2 & . & . & . & . & . & . & . & . \\ . & . & 6 & 16 & 2 & . & . & . & . & . & . & . \\ . & . & . & 20 & 25 & 2 & . & . & . & . & . & . \\ . & . & . & 8 & 50 & 36 & 2 & . & . & . & . & . \\ . & . & . & . & 35 & 105 & 49 & 2 & . & . & . & . \\ . & . & . & . & 10 & 112 & 196 & 64 & 2 & . & . & . \\ . & . & . & . & . & 54 & 294 & 336 & 81 & 2 & . & . \\ . & . & . & . & . & 12 & 210 & 672 & 540 & 100 & 2 & . \\ . & . & . & . & . & . & 77 & 660 & 1386 & 825 & 121 & 2 \end{bmatrix} $$ We have $$ C^{-1}_{r,c} = D_{r,c}\cdot {B_{2(r-c)} \over c } $$ where $B_k$ are the Bernoulli-numbers and the column- and rowindices $c,r$ begin at $1$.
Perhaps an even better display is this:
$$ E= \Tiny \begin{bmatrix}
1 & . & . & . & . & . & . & . & . & . & . & . & . & . \\
. & 1 & . & . & . & . & . & . & . & . & . & . & . & . \\
. & 3 & 1 & . & . & . & . & . & . & . & . & . & . & . \\
. & . & 6 & 1 & . & . & . & . & . & . & . & . & . & . \\
. & . & 5 & 10 & 1 & . & . & . & . & . & . & . & . & . \\
. & . & . & 15 & 15 & 1 & . & . & . & . & . & . & . & . \\
. & . & . & 7 & 35 & 21 & 1 & . & . & . & . & . & . & . \\
. & . & . & . & 28 & 70 & 28 & 1 & . & . & . & . & . & . \\
. & . & . & . & 9 & 84 & 126 & 36 & 1 & . & . & . & . & . \\
. & . & . & . & . & 45 & 210 & 210 & 45 & 1 & . & . & . & . \\
. & . & . & . & . & 11 & 165 & 462 & 330 & 55 & 1 & . & . & . \\
. & . & . & . & . & . & 66 & 495 & 924 & 495 & 66 & 1 & . & . \\
. & . & . & . & . & . & 13 & 286 & 1287 & 1716 & 715 & 78 & 1 & . \\
. & . & . & . & . & . & . & 91 & 1001 & 3003 & 3003 & 1001 & 91 & 1
\end{bmatrix}$$
where $ C^{-1}_{r,c} = E_{r,c} \cdot 2 B_{2(r-c)} / r $ and it seems $E_{r,c} = \binom{r}{2(r-c)} $ where $ c \le r \lt 2c$ so it seems
$$ C^{-1}_{r,c} = B_{2(r-c)} \cdot \binom{r}{2(r-c)}\cdot \frac 2{2c-r+1} \qquad \text{ where } c \le r \lt 2c \tag 1$$
[added]
If I use this description of the composition of $C^{-1}$ then I can symbolically invert this again and get (the top-left range here)
$$ C= \Tiny \begin{bmatrix}
1/2 & . & . & . & . & . \\
0 & 1 & . & . & . & . \\
0 & -3 \; B_2 & 3/2 & . & . & . \\
0 & 18 \; B_2^2 & -9 \; B_2 & 2 & . & . \\
0 & -180 \; B_2^3+15 \; B_4 \; B_2 & 90 \; B_2^2-15/2 \; B_4 & -20 \; B_2 & 5/2 & . \\
0 & 2700 \; B_2^4-495 \; B_4 \; B_2^2 & -1350 \; B_2^3+495/2 \; B_4 \; B_2 & 300 \; B_2^2-30 \; B_4 & -75/2 \; B_2 & 3
\end{bmatrix} \tag 2$$
and where the first row/column should be deleted since they are out of the pattern. For example, using the third and the fourth row in this matrix we get
$$ \begin{array} {rrlll} S(5,x) &=& -3 \; B_2 \; S(1,x)^2 &+ 3/2 \; S(2,x)^2 \\
S(7,x) &=& 18 \; B_2^2 \; S(1,x)^2 &- 9 \; B_2 \; S(2,x)^2 &+ 2 \; S(3,x)^2 \end{array} \tag 3$$
Too long for a comment on @Gottfried's answer. (Note: This was written before the answer was edited to include notes about the matrix $E$.)
In matrix $D$, the main diagonal consists of $2$s, and the next consists of squares: $$\begin{align} D_{r,r\phantom{-1}} &= 2 & r \geq 1\\ D_{r,r-1} &= (r-1)^2 = \binom {r-1}{1}\frac{r-1}{1} & r \geq 3 \end{align}$$
The third diagonal consists of consecutive 4-dimensional pyramidal numbers (OEIS #A002415); the next, consecutive 6-dimensional square numbers (OEIS #A040977); the next, consecutive 8-dimensional square numbers (OEIS #A053347).
If I've done my index arithmetic correctly (please double-check), we have $$\begin{align} D_{r,r-2} &= \binom{r-1}{3}\frac{r-2}{2} & r \geq 5 \\ D_{r,r-3} &= \binom{r-1}{5}\frac{r-3}{3} & r \geq 7 \\ D_{r,r-4} &= \binom{r-1}{7}\frac{r-4}{4} & r \geq 9 \end{align}$$ where row and column indices start at $1$. The two entries of the final visible column satisfy $$D_{r,r-5} = \binom{r-1}{9}\frac{r-5}{5} \qquad r \geq 11$$
So, for $r\neq c$, the non-zero entries of $D$ appear to be $$D_{r,c} = \binom{r-1}{2(r-c)-1}\frac{c}{r-c}$$
This is not an answer but does not fit in a comment-box.
Another way to describe the relation between the sums and the sums of squares of various sums-of-like-powers makes use of the Eulerian-numbers, which converts the expressions for the sums-of-like-powers into polynomials. The sums of like powers can be expressed as $$\small \begin{array} {rclllll} s_0(n) = & 1 \cdot \binom{n}{1} \\ s_1(n) = & & 1 \cdot \binom{n+1}{2} \\ s_2(n) = & & 1 \cdot \binom{n+1}{3}& +1 \cdot \binom{n+2}{3} \\ s_3(n) = & & 1 \cdot \binom{n+1}{4}& +4 \cdot \binom{n+2}{4} & +1 \cdot \binom{n+3}{4} \\ s_4(n) = & & 1 \cdot \binom{n+1}{5}& +11 \cdot \binom{n+2}{5} & +11 \cdot \binom{n+3}{5} & +1 \cdot \binom{n+4}{5} \\ s_5(n) = & & 1 \cdot \binom{n+1}{6}& +26 \cdot \binom{n+2}{6}& +66 \cdot \binom{n+3}{6} & +26 \cdot \binom{n+4}{6} & +1 \cdot \binom{n+5}{6} \\ \vdots & &\vdots \end{array} $$ This looks a bit more regular than the expressions by Bernoulli-polynomials and possibly the relations between the sums-expressions of odd orders $s_{2k+1}(n)$ can be written using these polynomials nicely and with more ease than with the Bernoulli-polynomials, don't know.
The first of your equality would then be expressed as $$ \begin{array} {rcl} \sum_{k=1}^n k^3 &=&\left(\sum_{k=1}^n k \right)^2 \\ s_3(n) &=& s_1(n)^2 \\ 1 \cdot \binom{n+1}{4} +4 \cdot \binom{n+2}{4} +1 \cdot \binom{n+3}{4} &=& \left( 1 \cdot \binom{n+1}{2}\right)^2 \\ \end{array}$$ and which must then be expanded to reproduce the equality (but I've not the time at the moment to do this in the required length...)
