How do I calculate these sum-of-sum expressions in terms of the generalized harmonic number? I know that $$\sum_{m=2}^k\sum_{n=1}^{m-1}(nm)^{-s}=\frac 12((H_k^s)^2-H_k^{(2s)})$$ and $H_k^s=\sum_{n=1}^kn^{-s}$
But, how would I go about finding identities in terms of the harmonic number like this for similar expressions with more sums? For example, how would I express
$$\sum_{q=3}^k\sum_{m=2}^{q-1}\sum_{n=1}^{m-1}(nmq)^{-s}$$ or $$\sum_{r=4}^k\sum_{q=3}^{r-1}\sum_{m=2}^{q-1}\sum_{n=1}^{m-1}(nmqr)^{-s}$$
and so on, in terms of harmonic numbers like we can show with two sums? (Check out the derivation of the two-sum expression: What is the proof for this sum of sum generalized harmonic number?)
 A: Your sums are better described as 
$$H_{k,r}^{(s)}=\sum_{1\le n_1<n_2<\cdots<n_r\le k}n_1^{-s}n_2^{-s}\cdots n_r^{-s}.\tag{1}$$
In fact, this is a special case of elementary symmetric polynomials:
$$e_r(x_1,\cdots,x_k)=\sum_{1\le i_1<i_2<\cdots<i_r\le k}x_{i_1}x_{i_2}\cdots x_{i_r}. \tag{2}$$
Here we have $H_{k,r}^{(s)}=e_r (1^{-s},2^{-s},\cdots,k^{-s})$. To relate this to the usual $H_k^{(s)}$s, we must first consider power-sum symmetric polynomials:
$$p_r(x_1,\cdots,x_k)=\sum_{i=1}^k x_i^r. \tag{3}$$
We have the following formula relating the two types of symmetric polynomials:
$$e_r=\frac{1}{r!}\sum_{\sigma\in S_r}{\rm sgn}(\sigma)p_1^{c_1(\sigma)}p_2^{c_2(\sigma)}\cdots p_r^{c_r(\sigma)}. \tag{4}$$
Here $S_r$ is the symmetric group (permutations of $r$ things), ${\rm sgn}(\sigma)$ denotes sign of a permutation, and the function $c_i(\sigma)$ denotes the number of cycles of length $i$ in $\sigma$'s disjoint cycle notation. This is often related to the Newton-Girard identities, which can be used to recursively obtain $e_r$s in terms of $p_i$s if this direct formula in $(4)$ is too daunting. The identity states
$$p_n-e_1p_{n-1}+\cdots+(-1)^{n-1} e_{n-1}p_1+(-1)^n e_n=0. \tag{5}$$
To see why $(5)$ is true, consider Vieta's formulas:
$$(x-x_1)(x-x_2)\cdots(x-x_n)=x^n-e_1x^{n-1}+\cdots+(-1)^{n-1}e_{n-1}x+(-1)^ne_n \tag{6}$$
Plug $x_1,x_2,\cdots,x_n$ into the equation $(6)$ and then add up the resulting $n$ equations to get $(5)$.
Anyway, in $(5)$ you can isolate $e_n$ to write it in terms of $p_1,\cdots,p_n$ and $e_1,\cdots,e_{n-1}$, which means you can keep iterating this in order to find $e_n$ purely in terms of $p_1,\cdots,p_n$.
Some threads on MSE that might be of further interest in understanding $(4)$ are:


*

*A conjecture relating Multiple Zeta Values and the Polya Enumeration Theorem

*Generalizing Newton's identities: Trace formula for Schur functors

*Symmetric and exterior power of representation
If we apply $(4)$ to $(1)$ and $(2)$ with $x_n=n^{-s}$, we get the explicit formula
$$H_{k,r}^{(s)}=\frac{1}{r!}\sum_{\sigma\in S_r}{\rm sgn}(\sigma)\left(H_k^{(s)}\right)^{c_1(\sigma)} \left(H_k^{(2s)}\right)^{c_2(\sigma)}\cdots \left(H_k^{(rs)}\right)^{c_r(\sigma)}. \tag{7}$$
One can rewrite this as a sum over integer partitions of $r$ by counting the number of permutations in $S_r$ with given cycle type - see equation $(11)$ in my answer as "blue" (third bullet point link above). Examples:
$$H_{k,3}^{(s)}=\frac{1}{6}\left[\left(H_k^{(s)}\right)^3-3H_k^{(s)} H_k^{(2s)}+2H_k^{(3s)}\right], \tag{8}$$
$$H_{k,4}^{(s)}=\frac{1}{24}\left[\left(H_k^{(s)}\right)^4-6\left(H_k^{(s)}\right)^2H_k^{(2s)}+3\left(H_k^{(2s)}\right)^2+8H_k^{(s)}H_k^{(3s)}-6H_k^{(4s)}\right]. \tag{9}$$
See this Wikipedia section for these examples (for symmetric polynomials in general).
