There is the statement about $\text{U}(N)$ model:
It is not possible to solve $\text{U}(N)$ model even in large $N$ limit
I do not understand this statement. If I go to large $N$ limit I know that only planar diagrams contribute to correlators. For $\text{O}(N)$ model, it seems that we can sum all planar diagrams (according to the statement that $\text{O}(N)$ model is exactly solvable in large $N$ limit). Where a probem appears in $\text{U}(N)$ model?
I will review the exact solution of the vector O($N$) model, and then discuss what goes wrong when you try to extend this to fields transforming in the adjoint representation. We have the (Euclidean) partition function $$ \mathcal{Z} = \int \mathcal{D}\phi_{\alpha} \, \delta\left(\phi_{\alpha}^2 - N \right) \exp\left[ - \frac{1}{2g^2} \int d^d x \left(\partial_{\mu} \phi_{\alpha} \right)^2 \right]. $$ Here I'm taking $\phi_{\alpha}^2 = \sum_{\alpha = 1}^N \phi_{\alpha} \phi_{\alpha}$ etc. Now, using the identity $\int d^x e^{i x p} = 2 \pi \delta(p)$, we can write $$ \mathcal{Z} = \int \mathcal{D}\phi_{\alpha} \mathcal{D}\lambda \exp\left[ - \frac{1}{2g^2} \int d^d x \left(\partial_{\mu} \phi_{\alpha} \right)^2 - \frac{1}{2} \int d^d x \, i \lambda \left( \phi_{\alpha}^2 - N \right) \right] $$ (I always absorb any overall constants into the path integral measure). Then one can integrate out the original scalar fields, $$ \mathcal{Z} = \int \mathcal{D}\lambda \exp\left[ - \frac{N}{2} \mathrm{Tr} \log \left( - \partial_{\mu}^2 + i \lambda \right) + \frac{N}{2g^2} \int d^d x \, i \lambda \right], $$ where I did some field redefinitions to move the $g$-dependence to the last term.
From here, we can evaluate this at large $N$ using the usual steepest descent method. If you would like a review of this method, I like Vadim Kaplunovsky's lecture notes in this pdf. As discussed there, this path integral is dominated by the maximum of the argument of the exponential, which has a width $1/\sqrt{N}$. Therefore, it is natural to change integration variables to $$ i \lambda = m^2 + \frac{i \beta}{\sqrt{N}}, $$ where $m^2$ is chosen to be the minimum of the action: $$ \int \frac{d^d p}{(2 \pi)^d} \frac{1}{p^2 + m^2} = \frac{1}{g^2}. $$ From here, we can expand the action in powers of $1/N$ using the new integration variable $\beta$, and we find that we're left with a Gaussian integral which can be performed exactly: $$ Z = e^{-\mathcal{S}[m^2]} \int \mathcal{D}\beta \exp\left[ - \frac{1}{2} \int \frac{d^d p}{(2\pi)^d} \Pi(p) |\beta(p)|^2 + O(1/\sqrt{N}) \right], $$ where $$ \Pi(p) = \int \frac{d^d k}{(2\pi)^d} \frac{1}{(k^2 + m^2)((k+p)^2 + m^2)}. $$ Then we can drop terms higher order in $1/N$, and everything has been solved. We only needed to perform two loop integrals: one to determine $m^2$, and one in the definition of $\Pi(p)$. If you want correlation functions, you just repeat this procedure with insertions of fields in the path integral.
So what goes wrong if we have fields transforming in the adjoint? Let's just try repeating the procedure with Polyakov's "principal chiral field" model: $$ \mathcal{Z} = \int \mathcal{D}g_{ab} \, \delta\left( \sum_c g_{ac}g^{\ast}_{cb} - \delta_{ab} N \right) \exp\left[ - \frac{1}{g^2} \int d^d x \left|\partial_{\mu} g_{ab} \right|^2 \right]. $$ Now we need a matrix of fields $\lambda_{ab}$ to decouple the delta function, which should really be thought of as a product of $N^2$ delta functions, one for each component of the matrix equation of its argument. This leads to $$ \mathcal{Z} = \int \mathcal{D}\lambda_{ab} \mathcal{D}g_{ab} \exp\left[ - \frac{1}{g^2} \int d^d x \left|\partial_{\mu} g_{ab} \right|^2 - \frac{1}{2} \int d^d x \, i \lambda_{ab} \left( g_{ac}g^{\ast}_{cb} - \delta_{ab} N \right) \right]. $$ Then, integrating out the $g$ fields, $$ \mathcal{Z} = \int \mathcal{D}\lambda_{ab} \exp\left[ - N \, \mathrm{Tr} \log \left( - \partial^2 \delta_{ab} + \lambda_{ab} \right) + \frac{N}{2g^2} \int d^d x \, \mathrm{Tr} \, i \lambda_{ab} \right]. $$ We now perform the saddle-point expansion as before, but I'll be careful to include the extra contributions for large-$N$. We end up with $$ Z = e^{-\mathcal{S}[m^2]} \int \mathcal{D}\beta_{ab} \exp\left[ - \frac{1}{2} \int \frac{d^d p}{(2\pi)^d} \Pi_{abcd}(p) \beta_{ab}(-p)\beta_{cd}(p) \right] \\ \times \, \exp\left[- \frac{1}{\sqrt{N}}\int \frac{d^d p_1 d^d p_2 d^d p_3}{(2 \pi)^{3d}} K^{(3)}_{abcdef}(p_1,p_2,p_3) \beta_{ab}(p_1)\beta_{cd}(p_2)\beta_{ef}(p_3) + \cdots \right]. $$ As you can probably convince yourself, the $n$-point vertex of the $\beta_{ab}$'s comes with a power of $N^{-(n-2)/2}$. However, unlike the vector case, diagrams will also involve summation over the indices $ab$, which will generate powers on $N$ in the numerator! So for example, if I consider a correction to the $\beta$ propagator due to the diagram (picture from Srednicki's textbook):
Each vertex contributes a factor of $1/\sqrt{N}$, but now we also need to sum over the indices associated with $K^{(3)}$ and the internal propagators; I'll leave it to you to show that this contributes a factor of $N$. And indeed, you can similarly show that all planar diagrams will contribute at $N=\infty$, so we have a long way to go before we can solve this theory. In contrast, for the vector O($N$) model, this diagram contributes a $1/N$ correction to the $N=\infty$ solution.
