I don't think that I have full solution to this problem (and I don't think it yet exists), however, for $N=3$ I think it is possible to make a sufficient comment.
First, one can note that each unitary matrix $U\in U(N)$ can be uniquely presented in the following form: \begin{align} U = \begin{bmatrix} e^{i\theta_1} & 0 & 0 & \cdots & 0\\ 0 & e^{i\theta_2} 0 & 0 &\cdots &0 \\ \vdots & \vdots & \vdots &\ddots & \vdots\\ 0 & 0 & 0 & \cdots & e^{i\theta_n} \end{bmatrix} \tilde{U} \begin{bmatrix} 1 & 0 & 0 & \cdots & 0\\ 0 & e^{i\mu_2} 0 & 0 &\cdots &0 \\ \vdots & \vdots & \vdots &\ddots & \vdots\\ 0 & 0 & 0 & \cdots & e^{i\mu_n} \end{bmatrix}, \end{align} where $\tilde{U}_{1j} = |U_{1j}|, \, \tilde{U}_{i1} = |U_{i1}|, \, i,j\in \{1,N\}$. Indeed, first column and first row of $U$ have $2N-1$ phases, which can be turned to identities by the $2N-1$ phases of $\theta,\, \mu$. One just need to solve a linear system to do the decomposition.
So now one may ask yourself, is it possible to reconstruct the unitary matrix $\tilde{U}$, knowing the $|U_{i,j}|$ and the fact that first row and first column are real, non-negative? For $N=3$ it looks possible.
So, for $N=3$ the unitary matrix $\tilde{U}$ is presented in the following form: \begin{equation} \tilde{U} = \begin{bmatrix} u_{11} & u_{12} & u_{13} \\ u_{21} & u_{22}e^{i\phi_{22}} & u_{23}e^{i\phi_{23}} \\ u_{31} & u_{32}e^{i\phi_{32}} & u_{33}e^{i\phi} \end{bmatrix}, \end{equation} where $\phi_{ij}$ are the phases we want to reconstruct from the knowledge of $u_{ij}$ in the formlua above. From the orthogonality relations for the columns of $\tilde{U}$ one writes the non-linear system of complex equations (it is important to note that it is non-linear): \begin{equation} \begin{cases} u_{11}u_{12} + u_{21}u_{22}e^{i\phi_{22}} + u_{31}u_{32}e^{i\phi_{32}} = 0,\\ u_{11}u_{13} + u_{21}u_{23}e^{i\phi_{23}} + u_{31}u_{33}e^{i\phi_{33}} = 0,\\ u_{12}u_{13} + u_{22}u_{23}e^{i(-\phi_{22} + \phi_{23})} + u_{32}u_{33}e^{i(-\phi_{32} + \phi_{33})} = 0 \end{cases} \end{equation} It is important that the system is non-linear, so the linearization argument for phases $\phi_{ij}$ doesn't work here, because for $\phi_{ij}\rightarrow 0$ matrix $\tilde{U}$ immediately appears to be non-unitary. In fact such system has a finite number of solutions (in fact no more than two). It is proved just by school methods of solving this system, so I show it here.
I will just consider the first equation from the system, which gives me two real equations for real and imaginary parts, respectively: \begin{equation} \begin{cases} u_{21}u_{22}\sin(\phi_{22}) + u_{31}u_{32}\sin(\phi_{32}) = 0, \\ u_{11}u_{12} + u_{21}u_{22}\cos(\phi_{22}) + u_{31}u_{32}\cos(\phi_{32}) = 0 \end{cases}. \end{equation} Solving the above system we obtain: \begin{align} \cos(\phi_{22}) &= \dfrac{u^2_{31}u^2_{32} - u^2_{21}u^2_{22}-u_{11}^2u^2_{12}} {2u_{11}u_{12}u_{21}u_{22}}, \\ \cos(\phi_{33}) &= \dfrac{u^2_{21}u^2_{22} - u^2_{31}u^2_{32}-u_{11}^2u^2_{12}} {2u_{11}u_{12}u_{31}u_{32}}. \end{align} Note these exact solutions imply that number of solutions of $\phi_{ij}$ is no more than finite and the reconstruction formulas are direct. One can also note that the non-linear complex system remains valid under complex conjugation, so if $\phi_{ij}$ is a solution, than $-\phi_{ij}$ is also a solution.
Discussion So what does this result say about our problem? The above considerations are taken from the work:
Auberson, G., Andre Martin, and G. Mennessier. "On the reconstruction of a unitary matrix from its moduli." Communications in mathematical physics 140.3 (1991): 523-542.
In the introduction it is written, that for $N=3$ the number of solutions of the aforementioned system is unique, up to complex conjugation (exactly what we spoke about in the end).
From this, the following result follows:
Let $U,V$ be unitary matrices from $U(3)$, which have the same moduli of their elements. It implies that there exist $(\theta_1, \dots, \theta_3), \, (\mu_2, \mu_3)$ such that: \begin{equation} U = \mathrm{diag}(e^{i\theta_1},e^{i\theta},e^{i\theta_3})V \mathrm{diag}(1,e^{i\mu_1},e^{i\mu_2}) \end{equation} or \begin{equation} U = \mathrm{diag}(e^{i\theta_1},e^{i\theta},e^{i\theta_3})V^* \mathrm{diag}(1,e^{i\mu_1},e^{i\mu_2}), \end{equation} where * denotes the complex conjugation.
It implies that dimension of the group is 5, when $N=3$.
Further discussions: For $N >3$ the result of Auberson, G., claims that there are cases, when the number of the solutions of the related non-linear system is infinite (or even a continuum), which implies that action of the group cannot be represented by $2N-1$ phases and finite number of nice "sandwich formulas".