I don't think that I have a full solution to this problem (and I don't think that 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 &\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 &\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, the first column and the 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 and 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 can write down the the non-linear system of complex equations:
\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 can be 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 exactly one, 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. Then 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".
