Denoting by PN (A, θ) = det I − Ae −iθ the characteristic polynomial on the unit circle in the complex plane of an N × N random unitary matrix A, we calculate the kth moment, defined with respect to an average over A ∈ U (N ), of the random variable corresponding to the 2βth moment of PN (A, θ) with respect to the uniform measure dθ 2π , for all k, β ∈ N . These moments of moments have played an important role in recent investigations of the extreme value statistics of characteristic polynomials and their connections with log-correlated Gaussian fields. Our approach is based on a new combinatorial representation of the moments using the theory of symmetric functions, and an analysis of a second representation in terms of multiple contour integrals. Our main result is that the moments of moments are polynomials in N of degree k 2 β 2 − k + 1. This resolves a conjecture of Fyodorov & Keating [23] concerning the scaling of the moments with N as N → ∞, for k, β ∈ N. Indeed, it goes further in that we give a method for computing these polynomials explicitly and obtain a general formula for the leading coefficient.