Denoting by P-N(A,theta)=det(I-Ae(-i theta)) the characteristic polynomial on the unit circle in the complex plane of an NxN random unitary matrix A, we calculate the kth moment, defined with respect to an average over A is an element of U(N), of the random variable corresponding to the 2 beta th moment of P-N(A,theta) with respect to the uniform measure d theta/2 pi, for all k,beta is an element of 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)beta(2)-k+1. This resolves a conjecture of Fyodorov and Keating (Philos Trans R Soc A 372(2007):20120503, 2014) concerning the scaling of the moments with N as N ->infinity, for k,beta is an element of 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.